{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "401_Adams Bashforth Example.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "I3Ze-AJLbVQR" }, "source": [ "# Adams Bashforth\n", "#### John S Butler \n", "john.s.butler@tudublin.ie \n", "[Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)\n", "\n", "\n", "\n", "The Adams Bashforth method is an explicit multistep method. This notebook illustrates the 2 step Adams Bashforth method for a linear initial value problem, given by\n", "\n", "\\begin{equation}y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2) \\end{equation}\n", "with the initial condition\n", "\\begin{equation}y(0)=1.\\end{equation}\n", "The video below walks through the notebook." ] }, { "cell_type": "code", "metadata": { "id": "rMDVFNpRbVQT", "outputId": "68da9c6d-130f-42d7-8b71-34011d3645ce" }, "source": [ "from IPython.display import HTML\n", "HTML('')" ], "execution_count": null, "outputs": [ { "output_type": "stream", "text": [ "/Users/johnbutler/anaconda3/lib/python3.7/site-packages/IPython/core/display.py:717: UserWarning: Consider using IPython.display.IFrame instead\n", " warnings.warn(\"Consider using IPython.display.IFrame instead\")\n" ], "name": "stderr" }, { "output_type": "execute_result", "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": { "tags": [] }, "execution_count": 1 } ] }, { "cell_type": "markdown", "metadata": { "id": "w2lET9kIbVQa" }, "source": [ "## Python Libraries" ] }, { "cell_type": "code", "metadata": { "id": "gZOAG1ijbVQb" }, "source": [ "import numpy as np\n", "import math \n", "import pandas as pd\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "HhsS6PbbbVQf" }, "source": [ "### Defining the function\n", "\\begin{equation} f(t,y)=t-y.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "tGq5h9O2bVQg" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "qQa7YsftbVQj" }, "source": [ "## Discrete Interval\n", "Defining the step size $h$ from the interval range $a\\leq t \\leq b$ and number of steps $N$ \n", "\\begin{equation}h=\\frac{b−a}{N}.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=t_0+ih,\\end{equation}\n", "where $t0=a.$\n", "\n", "Here the interval is $0≤t≤2$ and number of step 4 \n", "\\begin{equation}h=\\frac{2−0}{4}=0.5.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=0+i0.5,\\end{equation}\n", "for $i=0,1,⋯,4.$" ] }, { "cell_type": "code", "metadata": { "id": "DMV0kXM5bVQk", "outputId": "938b863c-d658-4ab3-cc3c-8aceed8ea94c" }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=4\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t,0*t,'o:',color='red')\n", "plt.xlim((0,2))\n", "plt.title('Illustration of discrete time points for h=%s'%(h))" ], "execution_count": null, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 1.0, 'Illustration of discrete time points for h=0.5')" ] }, "metadata": { "tags": [] }, "execution_count": 3 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAEICAYAAADiGKj0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAeeElEQVR4nO3dfZxmdV3/8ddbFjBQuVeRu6XAErpTJ6y0WpVbE9aSag0M+Wn8MNSfPzMDSUWQwn4WZmqGiiKQYKi5qEggQg8liYFIAVtZNnQ37l1AFow7P78/zhn22uG6Zq7lmt05M7yej8f1mPM953vO+Z7v9Z1r3nNuZlJVSJIkqbueNNsNkCRJ0tQMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2qZXkNUm+3lOuJHvMZpsGSfKRJO+Yhf2+PsltSdYk2W6I+jcl2bedfnuSj234Vm44SXZtj32T2W7LVJJcl2TRRtjPTyf59yT3JnnTDGzvhCRnzUTbpPnGwKYntN5AsYG2/8kk7xlxG+sESYCqOrqqThqtdevdjk2Bvwb2r6qnVNUP1mf9qvrzqnrdhmldf6O+v5PXr6rvt8f+yMy0cMOoqr2r6tJh6o7YR28DLq2qp1bVBx7nNmZckm2TfD7JfUm+l+T3p6h7QpKH2iA+8frJjdleaRgGNmkWJVkw221YD88AngxcN9sNgTnXd/PVbjzO8bCB378PAQ/SjNnDgL9LsvcU9c9tg/jEa8UGbJv0uBjYpCEkuTTJ63rKj571SuPUJLcnuSfJt5L8bJKjaH5YvK39rf38tv5NSf40ybeA+5IsSHJskhvbS0vXJ/mttu5zgI8Av9Ju4+52/jpn7pL8YZLlSVYnWZrkWT3LKsnRSW5IcleSDyXJgOPcPMn7k9zcvt7fzns2sKytdneSSwas/+r2jMYPkhw/admjl7uSPDnJWW29u5NcmeQZ7bJtk3yi3f9dSf6pnb8oyaq2724FPtHOf3mSa9rtXJ7k59v5ZwK7Aue3ffe2dv4vt/XuTvIfgy4d9ls/ycK2Pxf0jIv3tNtbk+T8JNslOTvJD9vjWtizzZ9JclH7Pi1L8rv99t2z7b9I8m/tuPpCkm17lh+S5tLn3W3d5/Qs670UfUKSzyT5VDu+rksyNsUxDnxvJrXvEuDFwAfbdZ+dZKt2P3e04+DPkjyprf+aJN9I872yGjhhwKFv1q+tw0qyJfBK4B1Vtaaqvg4sBV69PtuROqeqfPl6wr6Am4B92+nXAF/vWVbAHu30pcDrepY9Whc4ALgK2BoI8Bxgx3bZJ4H39NnnNcAuwE+0834HeBbNL1G/B9zXs4112jV5u8BLgDuB5wGbA38L/Muk4/hi275dgTuAAwf0x4nAN4GnAzsAlwMntcsWtttaMGDdvYA1wK+37fhr4OGe/j0BOKud/t/A+cAWwCbA84Gntcu+BJwLbANsCvxGO39Ru733ttv/ifaYbwde0G7niLZ/N5/8/rblnYAfAC9r+3q/trzDdOOjXx+042I58FPAVsD1wHeBfYEFwKeAT7R1twRWAke2y57Xvm97D9j3pcB/Az/brvvZnv57djtG9mv76G1tOzbrM65PAP6nPeZNgL8AvjnFMQ58bwa0sff74lPAF4Cntn31XeC1PeP4YeCN7fH/RJ/tTdfWLwJ3D3h9sa3zXOBHk7b7VuD8AcdwAnAPsJrmbOHrZ/tzyZevfi/PsEmje4jmB9TPAKmq71TVLdOs84GqWllVPwKoqn+sqpur6sdVdS5wA7DPkPs/DDi9qq6uqgeA42jOyC3sqXNKVd1dVd8Hvgb84hTbOrGqbq+qO4B3M/yZiUNpfmj+S9uOdwA/HlD3IWA7mkD8SFVdVVU/TLIjcBBwdFXdVVUPVdVlPev9GHhXVT3Q9t0fAn9fVVe02zkDeAD45QH7PRz4clV9ue3ri4BxmoDweH2iqm6sqnuAC4Abq+riqnoY+EeaAAHwcuCmqvpEVT1cVVfThLBDp9j2mVV1bVXdR9Ofv5vmgYffA75UVRdV1UPA+2gC7K8O2M7X22N+BDgT+IUp9tn3vZmuE3radVxV3VtVNwF/xbrj5+aq+tv2+H+0vm2tqpdX1dYDXi9vqz2FJoD1uofme7Sfz9D8krUDzXh6Z5JXTXe80sZmYJNGVFWXAB+kuW/mtiSnJXnaNKut7C0k+YOey3p305xV2X7IJjwL+F5Pe9bQnDXaqafOrT3T99P8UJt2W+30swbU7bfuo8fVhoxBDyacCVwInNNe+vzLNA817AKsrqq7Bqx3R1X9T095N+CPJ/qt7btdpmjzbsDvTKr/ImDHIY+xn9t6pn/UpzzR17sBL5i078OAZ06x7d5x8j2as2nb89j3/Mdt3Z3ob/L7/+QMvods0Hszne2BzXjs+Olt0zrjfgba2s8aYPL339OAe/tVrqrr21+WHqmqy4G/YeoQLc0KA5s0nPtoLhFNWOeHbFV9oKqeD+xNc7nqTyYWDdjeo/OT7AZ8FHgDsF1VbQ1cS3N5daptTLiZJgxMbG9LmjMk/z3NetNui+YS6s1DrnsLTViaaMcWbTseoz1z9u6q2ovmrNDLgT+g+YG+bZKtB+xjcl+sBE6edKZli6r69BT1z5xUf8uqOmXI/Y1iJXDZpH0/papeP8U6u/RM70pz9utOHvuep637eN7zdY5xivdmOne27Zs8fnrbNFJ/Jrkg6z7N2fu6oK32XWBBkj17Vv0Fhn84olj7vSd1hoFNGs41wG8n2SLN32Z77cSCJL+U5AXtWYj7aO7BmfizD7cB0/2JgC1pfkjc0W7vSJozbBNuA3ZOstmA9f8BODLJLybZHPhz4Ir2ktT6+jTwZ0l2SLI98E5g2L+LdR7w8iQvatt6IgM+Y5K8OMnPtZfRfkjzg/6R9lLyBcCHk2yTZNMkvz7FPj8KHN32f5JsmeQ3k0xc/prc/2cBByc5IMkm7Q32i5LsPGD7w7x/w/oi8Ow0D2Zs2r5+qfdhgT4OT7JXG35PBM5rLxV+BvjNJC9tx90f01wKvvxxtGudYxz03ky3kZ52nZzkqe0vIm9h+PEzrao6qNZ9mrP3dVBb5z7gc8CJ7Xh4IbCY5szhYyRZ3I61JNkHeBPNfXhSpxjYpOGcSvNnAm4DzgDO7ln2NJrgcBfNJaAf0NxTBPBxYK/2Etg/9dtwVV1Pc6/Pv7bb/zngGz1VLqE5O3Brkjv7rP9VmvubPktzluungCWP6yjhPTT3dH0L+DZwdTtvWlV1HXAMTYC8haY/Vg2o/kyagPdD4DvAZaz9wf5qmpDwnzQPFLx5in2O09x39MF2f8tpbm6f8Bc0AfTuJG+tqpU0P7zfThOQV9KcDR30WbjO+oPaMYyquhfYn+a9uZnm0t/EAxSDnEnzgMmtNH9S5U3ttpbR3I/3tzRntg4GDq6qBx9H0yYf41TvzXTeSPNLywrg6zRj4fTH0aZR/RHNPX230/wS8vp2fJLk15Ks6am7hGbc3Evz0MR723shpU5J1Uye8ZckzYQkl9I8FTqn/zuEpJnhGTZJkqSOM7BJkiR1nJdEJUmSOs4zbJIkSR03J/958vbbb18LFy6c7WZIkiRN66qrrrqzqnYYZRtzMrAtXLiQ8fHx2W6GJEnStJJ8b/paU/OSqCRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUsfNSGBLcmCSZUmWJzm2z/LNk5zbLr8iycJJy3dNsibJW2eiPZIkSfPJyIEtySbAh4CDgL2AVyXZa1K11wJ3VdUewKnAeyctPxW4YNS2SJIkzUczcYZtH2B5Va2oqgeBc4DFk+osBs5op88DXpokAEleAawArpuBtkiSJM07MxHYdgJW9pRXtfP61qmqh4F7gO2SbAn8KfDu6XaS5Kgk40nG77jjjhlotiRJ0twwE4EtfebVkHXeDZxaVWum20lVnVZVY1U1tsMOOzyOZkqSJM1NC2ZgG6uAXXrKOwM3D6izKskCYCtgNfAC4NAkfwlsDfw4yf9U1QdnoF2SJEnzwkwEtiuBPZPsDvw3sAT4/Ul1lgJHAP8KHApcUlUF/NpEhSQnAGsMa5IkSesaObBV1cNJ3gBcCGwCnF5V1yU5ERivqqXAx4EzkyynObO2ZNT9SpIkPVGkOdE1t4yNjdX4+PhsN0OSJGlaSa6qqrFRtuF/OpAkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkddyMBLYkByZZlmR5kmP7LN88ybnt8iuSLGzn75fkqiTfbr++ZCbaI0mSNJ+MHNiSbAJ8CDgI2At4VZK9JlV7LXBXVe0BnAq8t51/J3BwVf0ccARw5qjtkSRJmm9m4gzbPsDyqlpRVQ8C5wCLJ9VZDJzRTp8HvDRJqurfq+rmdv51wJOTbD4DbZIkSZo3ZiKw7QSs7Cmvauf1rVNVDwP3ANtNqvNK4N+r6oEZaJMkSdK8sWAGtpE+82p96iTZm+Yy6f4Dd5IcBRwFsOuuu65/KyVJkuaomTjDtgrYpae8M3DzoDpJFgBbAavb8s7A54E/qKobB+2kqk6rqrGqGtthhx1moNmSJElzw0wEtiuBPZPsnmQzYAmwdFKdpTQPFQAcClxSVZVka+BLwHFV9Y0ZaIskSdK8M3Jga+9JewNwIfAd4DNVdV2SE5Mc0lb7OLBdkuXAW4CJP/3xBmAP4B1JrmlfTx+1TZIkSfNJqibfbtZ9Y2NjNT4+PtvNkCRJmlaSq6pqbJRt+J8OJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnquBkJbEkOTLIsyfIkx/ZZvnmSc9vlVyRZ2LPsuHb+siQHDLXDq66ChQvh7LNnovmaj84+uxkjT3qSY0XTc7xoWI4VrY92vDwfnj/qphaMuoEkmwAfAvYDVgFXJllaVdf3VHstcFdV7ZFkCfBe4PeS7AUsAfYGngVcnOTZVfXItDv+3vfgqKOa6cMOG/UwNJ+cfXYzNu6/vyk7VjQVx4uG5VjR+pg8XkY0E2fY9gGWV9WKqnoQOAdYPKnOYuCMdvo84KVJ0s4/p6oeqKr/Apa32xvO/ffD298OixbBWWetnbdoEZx7blO+556m/LnPNeU772zK55/flG+9tSl/5StNeeXKpnzxxU15xYqmfNllTXnZsqZ8+eVN+dprm/KVVzbla65pytdc05SvvLIpX3ttU7788qa8bFlTvuyyprxiRVO++OKmvHJlU/7KV5ryrbc25fPPb8p33tmUP/e5pnzPPU353HOb8sQAOeuspvzQQ035k59syhM++lHYd9+15Q9/GA46aG35b/4GDjlkbfl974NXvnJt+ZRTYMmSteWTToLDD19bfuc74cgj15aPO27tBxzAW98KxxyztvzmNzevCccc09SZcNRRzTYmHHlks48Jhx/erDP5G+T+++H445u2v+99a+cfckhzjBMOOqjpgwn77tv00YRFi5o+hKZPHXtzf+wdf3z/8XL00es/9k46aW15yZKmjRMce3N/7A0aK8ccM/ufe469ptylsddvvIxgJgLbTsDKnvKqdl7fOlX1MHAPsN2Q6wKQ5Kgk40nG11mwcmW/6noim/hGmuz739+47dDcMGhcrFmzcduh7hs0VgZ95uiJbYZ/5qSqRttA8jvAAVX1urb8amCfqnpjT53r2jqr2vKNNGfSTgT+tarOaud/HPhyVX12qn2OJfVoatttN7jpppGOQfPMwoXNpYrJHCvqx/GiYTlWtD56xssYMF6VUTY3E2fYVgG79JR3Bm4eVCfJAmArYPWQ6w62xRZw8snr32LNbyef3IyNXo4VDeJ40bAcK1of/cbLCGYisF0J7Jlk9ySb0TxEsHRSnaXAEe30ocAl1ZzaWwosaZ8i3R3YE/i3ofa6225w2mne6KnHOuywZmzsthskjhVNzfGiYTlWtD56x8sMGPmSKECSlwHvBzYBTq+qk5OcCIxX1dIkTwbOBJ5Lc2ZtSVWtaNc9HvhfwMPAm6vqgun2NzY2VuPj49NVkyRJmnVJrqqqsZG2MROBbWMzsEmSpLliJgKb/+lAkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxIwW2JNsmuSjJDe3XbQbUO6Ktc0OSI9p5WyT5UpL/THJdklNGaYskSdJ8NeoZtmOBr1bVnsBX2/I6kmwLvAt4AbAP8K6eYPe+qvoZ4LnAC5McNGJ7JEmS5p1RA9ti4Ix2+gzgFX3qHABcVFWrq+ou4CLgwKq6v6q+BlBVDwJXAzuP2B5JkqR5Z9TA9oyqugWg/fr0PnV2Alb2lFe18x6VZGvgYJqzdJIkSeqxYLoKSS4Gntln0fFD7iN95lXP9hcAnwY+UFUrpmjHUcBRALvuuuuQu5YkSZr7pg1sVbXvoGVJbkuyY1XdkmRH4PY+1VYBi3rKOwOX9pRPA26oqvdP047T2rqMjY3VVHUlSZLmk1EviS4FjminjwC+0KfOhcD+SbZpHzbYv51HkvcAWwFvHrEdkiRJ89aoge0UYL8kNwD7tWWSjCX5GEBVrQZOAq5sXydW1eokO9NcVt0LuDrJNUleN2J7JEmS5p1Uzb2ri2NjYzU+Pj7bzZAkSZpWkquqamyUbfifDiRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6riRAluSbZNclOSG9us2A+od0da5IckRfZYvTXLtKG2RJEmar0Y9w3Ys8NWq2hP4alteR5JtgXcBLwD2Ad7VG+yS/DawZsR2SJIkzVujBrbFwBnt9BnAK/rUOQC4qKpWV9VdwEXAgQBJngK8BXjPiO2QJEmat0YNbM+oqlsA2q9P71NnJ2BlT3lVOw/gJOCvgPun21GSo5KMJxm/4447Rmu1JEnSHLJgugpJLgae2WfR8UPuI33mVZJfBPaoqv+bZOF0G6mq04DTAMbGxmrIfUuSJM150wa2qtp30LIktyXZsapuSbIjcHufaquART3lnYFLgV8Bnp/kprYdT09yaVUtQpIkSY8a9ZLoUmDiqc8jgC/0qXMhsH+SbdqHDfYHLqyqv6uqZ1XVQuBFwHcNa5IkSY81amA7BdgvyQ3Afm2ZJGNJPgZQVatp7lW7sn2d2M6TJEnSEFI1924HGxsbq/Hx8dluhiRJ0rSSXFVVY6Nsw/90IEmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOi5VNdttWG9J7gWWzXY7OmZ74M7ZbkQH2S/92S/92S+PZZ/0Z7/0Z7/099NV9dRRNrBgplqykS2rqrHZbkSXJBm3Tx7LfunPfunPfnks+6Q/+6U/+6W/JOOjbsNLopIkSR1nYJMkSeq4uRrYTpvtBnSQfdKf/dKf/dKf/fJY9kl/9kt/9kt/I/fLnHzoQJIk6Ylkrp5hkyRJesIwsEmSJHVcpwJbkgOTLEuyPMmxfZZvnuTcdvkVSRb2LDuunb8syQEbs90b2hD98pYk1yf5VpKvJtmtZ9kjSa5pX0s3bss3rCH65TVJ7ug5/tf1LDsiyQ3t64iN2/INZ4g+ObWnP76b5O6eZfN5rJye5PYk1w5YniQfaPvtW0me17Nsvo6V6frksLYvvpXk8iS/0LPspiTfbsfKyH+uoEuG6JdFSe7p+V55Z8+yKb//5rIh+uVPevrk2vbzZNt22bwcL0l2SfK1JN9Jcl2S/9Onzsx9tlRVJ17AJsCNwE8CmwH/Aew1qc4fAR9pp5cA57bTe7X1Nwd2b7ezyWwf00bslxcDW7TTr5/ol7a8ZraPYRb75TXAB/usuy2wov26TTu9zWwf08bok0n13wicPt/HSntsvw48D7h2wPKXARcAAX4ZuGI+j5Uh++RXJ44VOGiiT9ryTcD2s30Ms9Qvi4Av9pm/Xt9/c+01Xb9MqnswcMl8Hy/AjsDz2umnAt/t83Noxj5bunSGbR9geVWtqKoHgXOAxZPqLAbOaKfPA16aJO38c6rqgar6L2B5u735YNp+qaqvVdX9bfGbwM4buY2zYZjxMsgBwEVVtbqq7gIuAg7cQO3cmNa3T14FfHqjtGyWVdW/AKunqLIY+FQ1vglsnWRH5u9YmbZPqury9pjhifO5MsxYGWSUz6TOW89+eUJ8tlTVLVV1dTt9L/AdYKdJ1Wbss6VLgW0nYGVPeRWPPfBH61TVw8A9wHZDrjtXre+xvZYmzU94cpLxJN9M8ooN0cBZMmy/vLI9DX1ekl3Wc925Zujjai+b7w5c0jN7vo6VYQzqu/k6VtbX5M+VAv45yVVJjpqlNs2mX0nyH0kuSLJ3O8+xAiTZgiZ4fLZn9rwfL2lu0XoucMWkRTP22dKlf02VPvMm/82RQXWGWXeuGvrYkhwOjAG/0TN716q6OclPApck+XZV3bgB2rmxDdMv5wOfrqoHkhxNc3b2JUOuOxetz3EtAc6rqkd65s3XsTKMJ+Jny1CSvJgmsL2oZ/YL27HydOCiJP/ZnoF5Irga2K2q1iR5GfBPwJ44ViYcDHyjqnrPxs3r8ZLkKTQB9c1V9cPJi/us8rg+W7p0hm0VsEtPeWfg5kF1kiwAtqI5RTvMunPVUMeWZF/geOCQqnpgYn5V3dx+XQFcSvMbwHwwbb9U1Q96+uKjwPOHXXeOWp/jWsKkSxbzeKwMY1DfzdexMpQkPw98DFhcVT+YmN8zVm4HPs/8uQVlWlX1w6pa005/Gdg0yfY8wcdKj6k+W+bdeEmyKU1YO7uqPtenysx9tsz2TXs9N+YtoLnpbnfW3rC596Q6x7DuQwefaaf3Zt2HDlYwfx46GKZfnktzs+uek+ZvA2zeTm8P3MA8uQl2yH7ZsWf6t4BvttPbAv/V9s827fS2s31MG6NP2no/TXMTcJ4IY6XnGBcy+Eby32TdG4P/bT6PlSH7ZFea+4F/ddL8LYGn9kxfDhw428eyEfvlmRPfOzTB4/vtuBnq+28uv6bql3b5xEmULZ8I46V93z8FvH+KOjP22dKZS6JV9XCSNwAX0jxtc3pVXZfkRGC8qpYCHwfOTLKcZlAsade9LslngOuBh4Fjat1LPXPWkP3y/4CnAP/YPIPB96vqEOA5wN8n+THN2dRTqur6WTmQGTZkv7wpySE0Y2I1zVOjVNXqJCcBV7abO7HWPX0/Jw3ZJ9DcEHxOtZ8arXk7VgCSfJrm6b7tk6wC3gVsClBVHwG+TPM013LgfuDIdtm8HCswVJ+8k+Ye4Q+3nysPV9UY8Azg8+28BcA/VNVXNvoBbCBD9MuhwOuTPAz8CFjSfi/1/f6bhUPYIIboF2h+Mf7nqrqvZ9X5PF5eCLwa+HaSa9p5b6f5ZWfGP1v811SSJEkd16V72CRJktSHgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HH/H4txWz2MwbFSAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "tags": [], "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "wdFMNXnSbVQn" }, "source": [ "## Exact Solution\n", "THe initial value problem has the exact solution\n", "\\begin{equation}y(t)=2e^{-t}+t-1.\\end{equation}\n", "The figure below plots the exact solution." ] }, { "cell_type": "code", "metadata": { "id": "sGmWzAKrbVQn", "outputId": "f3999685-7c03-448c-869f-8b370ecb7dbe" }, "source": [ "IC=1 # Intial condtion\n", "y=(IC+1)*np.exp(-t)+t-1\n", "fig = plt.figure(figsize=(6,4))\n", "plt.plot(t,y,'o-',color='black')\n", "plt.title('Exact Solution ')\n", "plt.xlabel('time')" ], "execution_count": null, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 0, 'time')" ] }, "metadata": { "tags": [] }, "execution_count": 4 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXwU9f3H8dcH5BYEOSxGk4ilKgqipLZSEfAEpOCFgAEhIocVFRG8ULHSKFJu5TAoWGmERDyhCOXWQkEiCN6KchgBQZCjIOHI9/fHbvyFkJAEdnf2eD8fj324OzO7884wfvabz0xmzDmHiIhEvjJeBxARkcBQQRcRiRIq6CIiUUIFXUQkSqigi4hECRV0EZEooYIucgLMbIOZXXOC7403s/+ZWdlA55LYpoIuIeEvgL/4C1ne44Ugrq+FmWUXs8xZZvaGmf1kZrvN7BMz6x6ELEcVf+fcJufcqc65I4Fel8S2U7wOIDHlz865+V6HyGcqsAZIAHKAhsBvPE0kchI0QhfPmdkEM5uR7/VzZrbAfGqY2Swz225mP/ufn5Vv2dPNbIqZbfbPf9vMqgDvAWfm+23gzEJW/XvgFefcPufcYefcaufce/k+u52ZfWZmu8xssZldUET+V8zsb/le//rbgZlNBeKBmf4cD5lZopk5MzvFv8yZZvaume00s3Vm1jPfZz1lZplm9qqZ7fXnSTrRbS3RTQVdwsGDQCMz625mzYAeQDfnuy5FGWAKvlF0PPALkL9VMxWoDFwI1AFGOef2Aa2Bzf7WxqnOuc2FrHc5MM7MOplZfP4ZZvY7YBrQD6gNzMZXlMuX5gdzznUFNuH77eRU59ywQhabBmQDZwK3As+Y2dX55rcDpgPVgXcL/Pwiv1JBl1B62z/azXv0BHDO7Qe6ACOBfwL3Ouey/fN2OOfecM7td87tBVKB5gBmVhdf4e7jnPvZOXfIObekFHk6AB8ATwDrzexjM/u9f15H4F/OuXnOuUPAcKAS0PQkt8FRzOxs4ArgYefcAefcx8BLQNd8i/3HOTfb33OfClwcyAwSPVTQJZRudM5Vz/eYlDfDOfch8B1gQGbedDOrbGYvmtlGM9sDvA9U958hcjaw0zn384mE8X8JPOKcuxA4A/gY35eO4Rstb8y3bC7wPRB3Ius6jjPx/Qx7803bWGA9W/M93w9UzGvXiOSngi5hwczuASoAm4GH8s16EDgP+INzrhpwZd5b8BXY082seiEfWarLiDrnfsI3Cj8TON2fIyFfPsP3BfJDIW/fh6/tk6fggdXjZdmM72eomm9afBHrETkuFXTxnL9f/Td8bZeuwENm1tg/uyq+vvkuMzsdGJz3PufcFnwHP8f7D56WM7O8gv8jUNPMTjvOep8zs4vM7BR/Qb0bWOec24Hvt4QbzOxqMyuH74slB1hWyEd9DLTxH6D9Db6+e34/AvUKy+Cc+97/mc+aWUUza4TvGEJ6UblFiqKCLqGUd6ZH3uMtf+vgn8Bzzrk1zrlvgMeAqWZWARiNr3f9E76DmHMKfGZX4BDwJbANfzF1zn2J72Djd/5+fWFnuVQG3gJ24Wv3JOA7AIlz7it8XzDP+9f9Z3wHNg8W8jl5pz9uAP4NZBSY/yzwuD/HgELe3xlIxDdafwsY7JybV8hyIsdlusGFiEh00AhdRCRKqKCLiEQJFXQRkSihgi4iEiU8++OEWrVqucTERK9WLyISkT766KOfnHO1C5vnWUFPTEwkKyvLq9WLiEQkM9tY1Dy1XEREooQKuohIlFBBFxGJEiroIiJRQgVdRCRKqKCLiIRIeno6iYmJlClThsTERNLTA3tRTV0kX0QkBNLT0+nVqxf79+8HYOPGjfTq1QuA5OTkgKxDI3QRkRAYNGjQr8U8z/79+xk0aFDA1qGCLiISAps2bSrV9BOhgi4iEgJxcYXfjjY+Pj5g61BBFxEJgXPPPfeYaZUrVyY1NTVg61BBFxEJsjlz5rBkyRLatWtHQkICZkZCQgJpaWkBOyAKOstFRCSo9uzZQ8+ePbngggvIyMigYsWKQVuXCrqISBANHDiQzZs3s2zZsqAWc1DLRUQkaObPn09aWhoPPvggf/jDH4K+PnPOBX0lhUlKSnK6HrqIRKu9e/fSsGFDKlSowMcff0ylSpUC8rlm9pFzLqmweWq5iIgEwcMPP8ymTZv4z3/+E7BiXhy1XEREAmzRokVMmDCBfv360bRp05CtVy0XEZEA2rdvHw0bNqRs2bKsWbOGypUrB/Tz1XIREQmRRx99lA0bNrBkyZKAF/PiqOUiIhIgH3zwAc8//zx9+/alWbNmIV+/Wi4iIgGwf/9+Lr74Yo4cOcInn3xClSpVgrIetVxERILs8ccfZ926dSxcuDBoxbw4armIiJykZcuWMXr0aO6++25atmzpWQ4VdBGRk/DLL7+QkpJCfHw8zz33nKdZ1HIRETkJgwcP5uuvv2bevHlUrVrV0ywaoYuInKAVK1YwYsQIevbsyTXXXON1HBV0EZETceDAAVJSUoiLi2P48OFexwHUchEROSFPP/00X3zxBXPmzKFatWpexwE0QhcRKbWsrCyGDRvGnXfeyfXXX+91nF8VW9DNbLKZbTOzT4uYn2xma/2PZWZ2ceBjioiEh5ycHFJSUjjjjDMYMWKE13GOUpIR+itAq+PMXw80d841AoYAaQHIJSISllJTU/n0009JS0ujevXqXsc5SrE9dOfc+2aWeJz5y/K9XA6cdfKxRETCz+rVq3nmmWe44447uOGGG7yOc4xA99B7AO8VNdPMeplZlpllbd++PcCrFhEJnoMHD5KSkkLt2rUZNWqU13EKFbCzXMysJb6CfkVRyzjn0vC3ZJKSkry5KpiIyAkYOnQoa9as4e233+b000/3Ok6hAlLQzawR8BLQ2jm3IxCfKSISLtauXcuQIUO4/fbbad++vddxinTSLRcziwfeBLo6574++UgiIuHj0KFDpKSkcPrppzN27Fiv4xxXsSN0M5sGtABqmVk2MBgoB+Ccmwg8CdQExpsZwOGirtUrIhJphg0bxqpVq3jjjTeoWbOm13GOSze4EBEpwmeffcall17KjTfeSEZGhtdxgOPf4EJ/KSoiUojDhw+TkpJCtWrVeOGFF7yOUyK6louISCFGjBjBypUrycjIoHbt2l7HKRGN0EVECvjyyy8ZPHgwN998Mx06dPA6TompoIuI5HPkyBFSUlKoUqUK48ePx3+yR0RQy0VEJJ/Ro0ezfPly0tPTOeOMM7yOUyoaoYuI+H399dc8/vjjtGvXjs6dO3sdp9RU0EVE8LVa7rzzTipWrMjEiRMjqtWSRy0XERHghRdeYOnSpfzjH/+gbt26Xsc5IRqhi0jMW7duHY8++iht2rSha9euXsc5YSroIhLTcnNz6dGjB+XLlyctLS0iWy151HIRkZg2YcIE3n//fV5++WXi4uK8jnNSNEIXkZi1fv16Hn74Ya6//npSUlK8jnPSVNBFJCY557jrrrsoU6ZMxLda8qjlIiIxKS0tjYULF/Liiy8SHx/vdZyA0AhdRGLOxo0bGTBgANdccw09e/b0Ok7AqKCLSExxztGzZ0+cc0yaNCkqWi151HIRkZgyefJk5s2bx7hx40hMTPQ6TkBphC4iMSM7O5v+/fvTokUL+vTp43WcgFNBF5GY4JyjV69eHD58mJdffpkyZaKv/KnlIiIx4dVXX+W9995jzJgx1KtXz+s4QRF9X1EiIgVs3ryZfv360axZM/r27et1nKBRQReRqOaco3fv3hw4cCBqWy151HIRkaiWnp7OrFmzGDlyJPXr1/c6TlBF71eViMS8rVu3ct9993H55Zdz3333eR0n6FTQRSQqOee4++672b9/P5MnT6Zs2bJeRwo6tVxEJCplZGTw9ttvM2zYMM4//3yv44SERugiEnW2bdtG3759ueyyy+jfv7/XcUKm2IJuZpPNbJuZfVrE/PPN7L9mlmNmAwIfUUSkdO655x727t3LlClTYqLVkqckI/RXgFbHmb8TuA8YHohAIiInY8aMGcyYMYOnnnqKBg0aeB0npIot6M659/EV7aLmb3POrQQOBTKYiEhp/fTTT/zlL3+hSZMmDBw40Os4IRfSg6Jm1gvoBUTNBeVFJHzce++97Nq1iwULFnDKKbF3zkdID4o659Kcc0nOuaTatWuHctUiEuXefvttpk+fzhNPPEHDhg29juMJneUiIhFv586d9OnTh8aNG/PII494Hcczsfc7iYhEnfvvv58dO3YwZ84cypUr53UczxRb0M1sGtACqGVm2cBgoByAc26imf0GyAKqAblm1g9o4JzbE7TUIiJ+M2fO5J///CdPPvkkjRs39jqOp8w558mKk5KSXFZWlifrFpHo8PPPP3PhhRdSq1YtsrKyKF++vNeRgs7MPnLOJRU2Ty0XEYlY/fv3Z9u2bcyaNSsminlxdFBURCLSe++9xyuvvMLDDz/MpZde6nWcsKCWi4hEnN27d3PhhRdy2mmnsWrVKipUqOB1pJBRy0VEosqAAQPYsmULb731VkwV8+Ko5SIiEeXf//43L730EgMGDOD3v/+913HCilouIhIx9u7dy0UXXUTlypVZvXo1FStW9DpSyKnlIiJR4aGHHuL7779n6dKlMVnMi6OWi4hEhIULFzJx4kT69+/P5Zdf7nWcsKSWi4iEvf/97380bNiQcuXKsWbNGipVquR1JM+o5SIiEe3RRx9l48aNvP/++zFdzIujlouIhLUlS5bwwgsvcN9993HFFVd4HSesqeUiImFr3759XHzxxTjnWLt2LVWqVPE6kufUchGRiDRo0CC+/fZbFi1apGJeAmq5iEhYWrp0KWPHjuWee+6hRYsWXseJCCroIhJ2fvnlF1JSUkhISGDo0KFex4kYarmISNh54okn+Oabb1iwYAGnnnqq13EihkboIhJWli9fzqhRo+jduzdXXXWV13Eiigq6iISNAwcOkJKSQlxcHMOGDfM6TsRRy0VEwsZTTz3Fl19+ydy5c6lWrZrXcSKORugiEhZWrlzJ3//+d3r06MF1113ndZyIpIIuIp7Lycmhe/fu1K1blxEjRngdJ2Kp5SIinhsyZAiff/45//rXvzjttNO8jhOxNEIXEU+tWrWKoUOH0q1bN9q0aeN1nIimgi4injl48CDdu3enTp06jBo1yus4EU8tFxHxzDPPPMMnn3zCu+++S40aNbyOE/E0QhcRT6xZs4bU1FSSk5P585//7HWcqFBsQTezyWa2zcw+LWK+mdlYM1tnZmvN7NLAx/RJT08nMTGRMmXKkJiYSHp6erBWJSJBdOjQIbp3707NmjUZM2aM13GiRklG6K8ArY4zvzVQ3//oBUw4+VjHSk9Pp1evXmzcuBHnHBs3bqRXr14q6iIR6LnnnuPjjz9mwoQJ1KxZ0+s4UaNEN7gws0RglnPuokLmvQgsds5N87/+CmjhnNtyvM8s7Q0uEhMT2bhx4zHTExIS2LBhQ4k/R0S89cknn9CkSRNuueUWpk2b5nWciHO8G1wEooceB3yf73W2f1phQXqZWZaZZW3fvr1UK9m0aVOppotI+Dl8+DApKSlUr16d559/3us4UScQBd0KmVbosN85l+acS3LOJdWuXbtUK4mPjy/VdBEJP8OHD+ejjz5i3Lhx1KpVy+s4UScQBT0bODvf67OAzQH43KOkpqZSuXLlY6bff//9gV6ViATB559/zuDBg7n11lvp0KGD13GiUiAK+rvAHf6zXf4I7C6uf34ikpOTSUtLIyEhATMjLi6OSpUq8c9//pOcnJxAr05EAujIkSPceeedVK1alXHjxnkdJ2qV5LTFacB/gfPMLNvMephZHzPr419kNvAdsA6YBPwlWGGTk5PZsGEDubm5ZGdnM23aNFatWsWAAQOCtUoRCYBRo0axYsUKnn/+eerUqeN1nKhVorNcgqG0Z7kUpX///owaNYoZM2Zwyy23BCCZiATSV199xcUXX0zr1q158803MSvssJuU1PHOcon4gn7w4EGaNWvGl19+yerVq6lXr14A0olIIBw5coQrr7ySL774gs8//5zf/OY3XkeKeME+bdFT5cuXZ/r06ZQpU4aOHTuqny4SRsaOHcuyZcsYO3asinkIRHxBBzjnnHOYMmUKWVlZPPTQQ17HERFg3bp1DBo0iLZt25KcnOx1nJgQFQUd4MYbb+T+++9n7NixvPXWW17HEYlpubm53HnnnZQvX56JEyeqbx4iUVPQAYYNG0ZSUhIpKSmsX7/e6zgiMWvcuHF88MEHjB49mri4Qv9wXIIgqgp6+fLlycjIAKBjx44cPHjQ40Qisee7777jkUceoXXr1nTr1s3rODElqgo6QL169Zg8eTIrV67k4Ycf9jqOSEzJzc2lR48elC1blhdffFGtlhCLuoIOcPPNN3PvvfcyevRo3nnnHa/jiMSMF198kcWLFzNy5EjOPvvs4t8gARXx56EXJScnhz/96U98++23rF69msTExKCtS0Rgw4YNNGzYkMsvv5y5c+dqdB4kUX0eelEqVKhARkYGubm5dOrUSf10kSByztGzZ08AJk2apGLukagt6ADnnnsuL7/8MitWrODRRx/1Oo5I1Ml/W8j58+fToUMHEhISvI4Vs6K6oAPceuut3HPPPYwcOZJ3333X6zgiUSP/bSHzZGRk6LaQHoraHnp+Bw4coGnTpmzYsIHVq1drBCESALotpDdisoeeX8WKFcnMzOTw4cN06tSJQ4cOeR1JJOLptpDhJyYKOsBvf/tbXnrpJZYvX85jjz3mdRyRiHb48GEqVapU6DzdFtI7MVPQAW677Tbuvvtuhg8fzqxZs7yOIxKRDh06ROfOndm/fz/lypU7al7lypVJTU31KJnEVEEHGDlyJI0bN6Zbt258//33XscRiSgHDx6kY8eOzJgxg5EjRzJlypRfbwuZkJBAWlqarqzooZg4KFrQN998w6WXXkqjRo1YvHjxMaMMETlWTk4OHTp0YObMmYwdO5Z7773X60gxKeYPihZUv359Jk2axLJly3j88ce9jiMS9g4cOMBNN93EzJkzmTBhgop5mIrJgg7QqVMnevfuzbBhw5g9e7bXcUTC1v79+2nXrh1z5sxh0qRJ9OnTp/g3iSditqCD707kjRo14o477iA7O9vrOCJhZ9++fbRt25b58+czZcoU7rrrLq8jyXHEdEGvVKkSmZmZ5OTk0KlTJw4fPux1JJGwsXfvXtq0acOSJUuYOnWqrm0eAWK6oAOcd955vPjiiyxdupQnnnjC6zgiYWHPnj20atWKpUuX8tprr+nMlQgR8wUd4Pbbb6dnz54MHTqUOXPmeB1HxFO7du3iuuuu48MPPyQjI4OOHTt6HUlKSAXdb8yYMTRs2JCuXbvyww8/eB1HxBM7d+7k2muvZdWqVcyYMYNbbrnF60hSCirofnn99F9++YXOnTurny4xZ8eOHVx99dWsXbuWN998k/bt23sdSUqpRAXdzFqZ2Vdmts7MHilkfoKZLTCztWa22MzOCnzU4Dv//POZOHEiH3zwAYMHD/Y6jkjIbN++nZYtW/LFF1/wzjvv0LZtW68jyQkotqCbWVlgHNAaaAB0NrMGBRYbDrzqnGsEPA08G+igodKlSxd69OjBs88+y9y5c72OIxJ0P/74Iy1btmTdunXMmjWLVq1aeR1JTlBJRuiXAeucc9855w4C04GCv4s1ABb4ny8qZH5EGTt2LBdeeCFdu3Zl8+bNXscRCZotW7bQokUL1q9fz+zZs7nmmmu8jiQnoSQFPQ7IfxWrbP+0/NYAeUdPbgKqmlnNk4/njcqVK5OZmcm+ffu4/fbb1U+XqJSdnU3z5s3Jzs5mzpw5tGjRwutIcpJKUtALu9trwSt6DQCam9lqoDnwA3BMFTSzXmaWZWZZ27dvL3XYULrggguYMGECS5Ys4a9//avXcUQCatOmTTRv3pytW7cyd+5cmjVr5nUkCYCSFPRs4Ox8r88CjupDOOc2O+duds5dAgzyT9td8IOcc2nOuSTnXFLt2rVPInZo3HHHHaSkpJCamsq8efO8jiMSEBs2bKB58+bs2LGD+fPn07RpU68jSYCUpKCvBOqb2TlmVh7oBBx1t2Uzq2VmeZ/1KDA5sDG98/zzz3PBBReQnJzMli1bvI4jclK+/fZbrrzySnbv3s2CBQu47LLLvI4kAVRsQXfOHQb6AnOBL4BM59xnZva0mbXzL9YC+MrMvgbOAKLmliVVqlTh9ddf/7WffuTIEa8jiZyQr7/+mubNm7N//34WLlxIkyZNvI4kARaTN7g4Ea+88gopKSk8+eST6qlLxPnyyy+56qqrOHz4MAsWLKBhw4ZeR5ITpBtcBED37t3p1q0bQ4YMYcGCBcW/QSRMfPbZZzRv3pzc3FwWL16sYh7FVNBLYdy4cZx//vkkJyezdetWr+OIFGvt2rW0aNGCsmXLsnjxYho0KPg3gRJNVNBLoUqVKmRmZrJnzx6Sk5PVT5ewtnr1alq2bEnFihVZsmQJ559/vteRJMhU0Evpoosu4oUXXmDhwoX87W9/8zqOSKGysrK46qqrOPXUU1myZAn169f3OpKEgAr6CUhJSaFr16789a9/ZeHChV7HETnKihUruOaaa6hRowZLliyhXr16XkeSEFFBPwFmxvjx4znvvPNITk7mxx9/9DqSCABLly7l2muvpVatWixZsoTExESvI0kIqaCfoFNPPZXMzEx27dpFly5d1E8Xz73//vtcf/311K1blyVLlnD22WcX/yaJKiroJ6Fhw4Y8//zzzJ8/n2eeecbrOBLDFi5cSOvWrYmPj2fx4sXExRW8fp7EAhX0k9SjRw+Sk5N56qmnWLx4sddxJAb9+9//5oYbbqBevXosWrSIunXreh1JPKKCfpLMjAkTJvDb3/6W22+/nW3btnkdSWLI7NmzadeuHb/73e9YuHAhZ5xxhteRxEMq6AFQtWpVXn/9dX7++We6dOlCbm6u15EkBsycOZObbrqJBg0asHDhQiLhCqYSXCroAdKoUSPGjBnDvHnzePbZiL0Dn0SIt956i5tvvpmLL76YBQsWULNmxN5PRgJIBT2AevbsSefOnXnyySd5//33vY4jUer111+nQ4cOJCUlMW/ePGrUqOF1JAkTKugBZGa8+OKLnHvuuXTu3JlwvyuTRJ5p06bRuXNn/vjHPzJ37lxOO+00ryNJGFFBD7CqVauSmZnJjh076Nq1q/rpEjBTp06lS5cuXHHFFcyZM4dq1ap5HUnCjAp6EDRu3JjRo0czd+5cnnvuOa/jSBSYPHky3bp1o2XLlsyePZtTTz3V60gShlTQg6R379507NiRJ554gg8++MDrOBLB0tLS6NGjB9dddx0zZ86kcuXKXkeSMKWCHiRmRlpaGueccw6dO3fmp59+8jqSRKBx48bRu3dvbrjhBt5++20qVarkdSQJYyroQVStWjUyMzPZvn07d9xxh/rpUipjxoyhb9++tG/fnjfeeIOKFSt6HUnCnAp6kF1yySWMGjWK9957j7///e9ex5EIMXz4cPr168fNN99MZmYmFSpU8DqSRAAV9BC4++676dChA4MGDWLp0qVex5Ew9+yzzzJw4EBuu+02pk+fTvny5b2OJBFCBT0EzIxJkyaRkJBAp06d2LFjh9eRJEw9/fTTPPbYY9x+++2kp6dTrlw5ryNJBFFBD5HTTjuNzMxMtm3bRrdu3dRPl6M453jyyScZPHgw3bp149VXX+WUU07xOpZEGBX0EGrSpAkjRozgX//6FyNGjPA6joQJ5xyPPfYYQ4YMoUePHkyePJmyZct6HUsikAp6iN1zzz3ccsstPProoyxbtszrOOIx5xwDBw5k6NCh9OnTh7S0NMqU0f+WcmK054SYmfHyyy8THx+vfnqMc87xwAMPMGLECPr27cv48eNVzOWkaO/xQF4/fevWrXTv3h3nnNeRJMRyc3Pp27cvY8aM4YEHHmDs2LGYmdexJMKVqKCbWSsz+8rM1pnZI4XMjzezRWa22szWmlmbwEeNLklJSQwfPpxZs2YxcuRIr+NICOXm5tKnTx/Gjx/PQw89xIgRI1TMJTCcc8d9AGWBb4F6QHlgDdCgwDJpwN3+5w2ADcV9bpMmTVysy83NdTfddJM75ZRT3H//+1+v40gIHD582KWkpDjADRo0yOXm5nodSSIMkOWKqKslGaFfBqxzzn3nnDsITAfaF/xeAPKu5XkasPkkvmNihpkxefJkzjrrLDp27MjOnTu9jiRBdOTIEVJSUpgyZQpPPfUUQ4YM0chcAqokBT0O+D7f62z/tPyeArqYWTYwG7i3sA8ys15mlmVmWbr5g0/16tXJyMhgy5YtpKSkqJ8epQ4fPkyXLl2YOnUqf/vb3xg8eLCKuQRcSQp6YXtdwarTGXjFOXcW0AaYambHfLZzLs05l+ScS9INbf/fZZddxrBhw3j33XcZPXq013EkwA4dOkTnzp2ZPn06zz33HIMGDfI6kkSpkhT0bODsfK/P4tiWSg8gE8A591+gIlArEAFjxf3330/79u15+OGH+fDDD72OIwFy8OBBOnbsyIwZMxg5ciQPPfSQ15EkipWkoK8E6pvZOWZWHugEvFtgmU3A1QBmdgG+gq6eSimYGVOmTOHMM8/ktttu4+eff/Y6kpyknJwcbr31Vt566y3Gjh3LAw884HUkiXLFFnTn3GGgLzAX+ALIdM59ZmZPm1k7/2IPAj3NbA0wDeju1AwutRo1apCRkcEPP/zAnXfeqX56BDtw4AA33XQTM2fOZMKECdx7b6GHlUQCq6jTX4L90GmLRRs5cqQD3OjRo72OIidg37597tprr3Vm5iZNmuR1HIkynORpixJi/fr1o127dgwcOJCVK1d6HUdKYd++fbRt25b58+czZcoU7rrrLq8jSQxRQQ9Def30unXrctttt7Fr1y6vI0kJ7N27lzZt2rBkyRKmTp1Kt27dvI4kMUYFPUydfvrpZGRkkJ2dTY8ePdRPD3N79uyhVatWLF26lNdee43k5GSvI0kMUkEPY3/84x8ZOnQob775Ji+88ILXcaQIu3bt4rrrruPDD4YdxjUAAAzBSURBVD8kIyODjh07eh1JYpQKepjr378/bdu2ZcCAAWRlZXkdRwrYuXMn1157LatWrWLGjBnccsstXkeSGKaCHubMjFdeeYUzzjiDjh07snv3bq8jid+OHTu4+uqrWbt2LW+++Sbt2xe8xJFIaKmgR4CaNWsyffp0Nm7cyF133aV+ehjYvn07LVu25IsvvuCdd96hbdu2XkcSUUGPFE2bNuXZZ59lxowZjB8/3us4MW3r1q20aNGCdevWMWvWLFq1auV1JBFABT2iPPjgg7Rp04b+/fuzatUqr+PEpM2bN9OiRQs2bNjA7Nmzueaaa7yOJPIrFfQIUqZMGf7xj39Qp04dbrvtNvbs2eN1pJiSnZ1NixYt+OGHH5gzZw4tWrTwOpLIUVTQI0ytWrWYPn06GzZsoGfPnuqnh8jGjRtp3rw5W7duZe7cuTRr1szrSCLHUEGPQH/6059ITU0lMzOTiRMneh0n6q1fv57mzZuzY8cO5s+fT9OmTb2OJFIoFfQINXDgQFq3bk2/fv1YvXq113Gi1rp162jevDl79uxhwYIFXHbZZV5HEimSCnqEyuun165dW/30IPn6669p3rw5+/fvZ+HChTRp0sTrSCLHpYIewWrXrs20adNYv349vXr1Uj89ANLT00lMTKRMmTJccMEF7N27l0WLFtG4cWOvo4kUSwU9wjVr1owhQ4aQkZFBWlqa13EiWnp6Or169WLjxo0458jNzeXQoUOsXbvW62giJWJejeqSkpKcrk0SGLm5ubRp04bFixezfPlyjSaLceDAAX744YdjHmlpaezfv/+Y5RMSEtiwYUPog4oUwsw+cs4lFTpPBT06bNu2jUsuuYQqVarw0UcfUbVqVa8jhZxzjh07dvxaoLOzswst3Dt37jzmvVWqVGHfvn2Ffq6ZkZubG+z4IiVyvIJ+SqjDSHDUqVOHadOm0bJlS3r37k16ejpm5nWsgMnJyWHz5s2FFui8x+bNm8nJyTnqfWZGnTp1iIuLIzExkSuuuIK4uLhjHtWqVeOcc85h48aNx6w7Pj4+VD+myElRQY8iV155JU8//TSPP/44LVu2pGfPnl5HKpZzjp9//rnIIp03yv7pp5+OeW+lSpV+LciXX355oYW6bt26lCtXrkRZUlNT6dWr11Ftl8qVK5Oamhqwn1ckmNRyiTK5ubm0atWKDz74gBUrVtCoUSPPshw6dIgtW7YU2frIexw4cOCY99auXbvQAh0XF8dZZ51FXFwc1atXD/hvIenp6QwaNIhNmzYRHx9Pamqq7j4kYUU99Bizbds2GjduTNWqVcnKygp4P905x+7du49bpH/44Qe2bdt2zKmUFSpU4MwzzzymOBccVVeoUCGgmUWihXroMaZOnTq89tprXH311bRq1Yrs7Gy+//77Eo04Dx8+zNatW4tsfeQ9CjsbpGbNmr8W5UsvvbTQ0XXNmjWjqrcvEk40Qo9it956K2+88cZR0ypWrEi/fv04//zzCx1V//jjj8ec0VGuXLmjRtWFPc4880wqVaoUyh9PJCZphB6jVq5cecy0AwcOMHTo0F9fV69e/dei3LBhw0LbILVq1aJMGf0Nmki4U0GPYt9//32h082Mr776iri4OCpXrhziVCISLCUadplZKzP7yszWmdkjhcwfZWYf+x9fm9muwEeV0irq/On4+Hjq16+vYi4SZYot6GZWFhgHtAYaAJ3NrEH+ZZxzDzjnGjvnGgPPA28GI6yUTmpq6jFFW+dVi0SvkozQLwPWOee+c84dBKYD7Y+zfGdgWiDCyclJTk4mLS2NhIQEzIyEhATS0tJ0XrVIlCpJDz0OyN+MzQb+UNiCZpYAnAMsPPloEgjJyckq4CIxoiQj9MJOGi7qXMdOwAzn3JFCP8isl5llmVnW9u3bS5pRRERKoCQFPRs4O9/rs4DNRSzbieO0W5xzac65JOdcUu3atUueUkREilWSgr4SqG9m55hZeXxF+92CC5nZeUAN4L+BjSgiIiVRbEF3zh0G+gJzgS+ATOfcZ2b2tJm1y7doZ2C6033QREQ8UaI/LHLOzQZmF5j2ZIHXTwUuloiIlJZn13Ixs+3AsXcTKJlawLEXyPZeuOaC8M2mXKWjXKUTjbkSnHOFHoT0rKCfDDPLKuriNF4K11wQvtmUq3SUq3RiLZeuuCQiEiVU0EVEokSkFvQ0rwMUIVxzQfhmU67SUa7SialcEdlDFxGRY0XqCF1ERApQQRcRiRJhV9BLcDONCmaW4Z+/wswS88171D/9KzO7PsS5+pvZ52a21swW+K88mTfvSL4bgBxz2YQg5+puZtvzrf+ufPO6mdk3/ke3EOcq8qYoQd5ek81sm5l9WsR8M7Ox/txrzezSfPOCub2Ky5Xsz7PWzJaZ2cX55m0ws0/82yugN+otQa4WZrY737/Xk/nmHXcfCHKugfkyferfp073zwvK9jKzs81skZl9YWafmdn9hSwT3P3LORc2D6As8C1QDygPrAEaFFjmL8BE//NOQIb/eQP/8hXwXcL3W6BsCHO1BCr7n9+dl8v/+n8ebq/uwAuFvPd04Dv/f2v4n9cIVa4Cy98LTA729vJ/9pXApcCnRcxvA7yH7yqjfwRWBHt7lTBX07z14bvZzIp88zYAtTzaXi2AWSe7DwQ6V4Fl/wwsDPb2AuoCl/qfVwW+LuT/x6DuX+E2Qi/JzTTaA//wP58BXG1m5p8+3TmX45xbD6zzf15IcjnnFjnn9vtfLsd3VcpgK+3NR/K7HpjnnNvpnPsZmAe08ihXyG6K4px7H9h5nEXaA686n+VAdTOrS3C3V7G5nHPL/OuF0O1fJdleRTmZfTPQuUKyfznntjjnVvmf78V37au4AosFdf8Kt4Je2M00Cm6QX5dxvguH7QZqlvC9wcyVXw9838J5KprvOvDLzezGAGUqTa5b/L/ezTCzvEshh8X2ssJvihKs7VUSRWUP5vYqrYL7lwP+bWYfmVkvD/JcbmZrzOw9M7vQPy0stpeZVcZXGN/INzno28t8reBLgBUFZgV1/yrRxblCqCQ30yhqmdLciKO0SvzZZtYFSAKa55sc75zbbGb1gIVm9olz7tsQ5ZoJTHPO5ZhZH3y/3VxVwvcGM1eewm6KEqztVRJe7F8lZmYt8RX0K/JN/pN/e9UB5pnZl/4RbCiswndtkf+ZWRvgbaA+YbK98LVbljrn8o/mg7q9zOxUfF8g/ZxzewrOLuQtAdu/wm2EXpKbafy6jJmdApyG71ev0tyIIxi5MLNrgEFAO+dcTt5059xm/3+/Axbj++YOSS7n3I58WSYBTUr63mDmyueYm6IEcXuVRFHZg7m9SsTMGgEvAe2dczvypufbXtuAtwhcq7FYzrk9zrn/+Z/PBsqZWS3CYHv5HW//Cvj2MrNy+Ip5unPuzUIWCe7+FegDAyd5UOEUfAcDzuH/D6RcWGCZezj6oGim//mFHH1Q9DsCd1C0JLkuwXcQqH6B6TWACv7ntYBvCNDBoRLmqpvv+U3Acvf/B2HW+/PV8D8/PVS5/Mudh+8AlYVie+VbRyJFH+S7gaMPWn0Y7O1Vwlzx+I4LNS0wvQpQNd/zZUCrEOb6Td6/H77CuMm/7Uq0DwQrl39+3mCvSii2l//nfhUYfZxlgrp/BWzjBvAfqQ2+o8PfAoP8057GN+oFqAi87t+5PwTq5XvvIP/7vgJahzjXfOBH4GP/413/9KbAJ/4d+hOgR4hzPQt85l//IuD8fO+9078d1wEpoczlf/0UMLTA+4K9vaYBW4BD+EZFPYA+QB//fAPG+XN/AiSFaHsVl+sl4Od8+1eWf3o9/7Za4/93HhTiXH3z7V/LyfeFU9g+EKpc/mW64ztRIv/7gra98LXBHLA2379Tm1DuX/rTfxGRKBFuPXQRETlBKugiIlFCBV1EJEqooIuIRAkVdBGRKKGCLjHDzKqb2V/8z880sxleZxIJJJ22KDHDf32NWc65izyOIhIU4XYtF5FgGgqca2Yf4/sL1AuccxeZWXfgRnyXfL0IGIHvrxu7AjlAG+fcTjM7F98fhdQG9gM9nXNfhv7HECmcWi4SSx4BvnXONQYGFph3EXA7vj9fTwX2O+cuAf4L3OFfJg241znXBBgAjA9JapES0ghdxGeR813Deq+Z7cZ3lUrw/Xl2I/8V9JoCr/suvw/4rhskEjZU0EV8cvI9z833Ohff/ydlgF3+0b1IWFLLRWLJXny3Bis157uu9Xoz6wC/3hvy4mLeJhJSKugSM5zvGuJL/TcW/vsJfEQy0MPM8q7UF7BbqokEgk5bFBGJEhqhi4hECRV0EZEooYIuIhIlVNBFRKKECrqISJRQQRcRiRIq6CIiUeL/AAHX1Are7O2LAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "tags": [], "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "V91-l3YzbVQq" }, "source": [ "## 2-step Adams Bashforth\n", "\n", "The general 2-step Adams Bashforth difference equation is\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{2}(3f(t_i,w_i)-f(t_{i-1},w_{i-1})). \\end{equation}\n", "For the specific intial value problem the 2-step Adams Bashforth difference equation is\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{2}(3(t_i-w_i)-(t_{i-1}-w_{i-1})). \\end{equation}\n", "\n", "for $i=0$ the difference equation is:\n", "\\begin{equation}w_{1} = w_{0} + \\frac{h}{2}(3(t_0-w_0)-(t_{-1}-w_{-1})), \\end{equation}\n", "this is not solvable as $w_{-1}$ is unknown.\n", "for $i=1$ the difference equation is:\n", "\\begin{equation}w_{2} = w_{1} + \\frac{h}{2}(3(t_1-w_1)-(t_{0}-w_{0})), \\end{equation}\n", "this is not solvable as $w_{1}$ is unknown, but it can be approximated using a one step method.\n", "Here, as the exact solution is known,\n", "\\begin{equation}w_1=2e^{-t_1}+t_1-1.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "A3mFbuf6bVQq" }, "source": [ "### Initial conditions\n", "w=np.zeros(len(t))\n", "w[0]=IC\n", "w[1]=y[1] # NEED FOR THE METHOD" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "XelJoDUcbVQs" }, "source": [ "### Loop" ] }, { "cell_type": "code", "metadata": { "id": "uX_cZAj_bVQt" }, "source": [ "for k in range (1,N):\n", " w[k+1]=w[k]+h/2.0*(3*myfun_ty(t[k],w[k])-myfun_ty(t[k-1],w[k-1])) " ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "apgJ8EZHbVQu" }, "source": [ "### Plotting solution" ] }, { "cell_type": "code", "metadata": { "id": "41bxcOCObVQu" }, "source": [ "def plotting(t,w,y):\n", " fig = plt.figure(figsize=(10,4))\n", " plt.plot(t,y, 'o-',color='black',label='Exact')\n", " plt.plot(t,w,'s:',color='blue',label='Adams-Bashforth')\n", " plt.xlabel('time')\n", " plt.legend()\n", " plt.show " ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "LDg9CIYhbVQw" }, "source": [ "The plot below shows the exact solution (black) and the 2 step Adams-Bashforth approximation (red) of the intial value problem" ] }, { "cell_type": "code", "metadata": { "id": "yRUWg92dbVQw", "outputId": "2a7378ef-e0fb-476b-c577-fd132e3803e2" }, "source": [ "plotting(t,w,y)" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEGCAYAAABB6hAxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3xO5//H8deViC12USP2JkYiNqlNxagWDWqmNVs/tUrt2Sqq+BYtWlK7SEtp7U2CUJuYqb2SGCHj+v1xxAxCTnLuJJ/n43E/cl/3ObnO54jxdp3rXEdprRFCCCGEEG/HzuoChBBCCCESMglTQgghhBCxIGFKCCGEECIWJEwJIYQQQsSChCkhhBBCiFhIZtWBs2TJovPmzWvV4YUQQgghYmzv3r3XtdZZo9tmWZjKmzcvfn5+Vh1eCCGEECLGlFLnXrZNLvMJIYQQQsSChCkhhBBCiFiQMCWEEEIIEQuWzZmKTlhYGIGBgYSGhlpdiohDKVOmJFeuXDg4OFhdihBCCBFrNhWmAgMDSZcuHXnz5kUpZXU5Ig5orblx4waBgYHky5fP6nKEEEKIWLOpy3yhoaFkzpxZglQippQic+bMMvoohBAiVrJnB6VefGXPHv+12FSYAiRIJQHyMxZCCBFbV6682edxyebClBBCCCFEQiJh6jn29vaUKVPm8WvcuHGm9e3v78/q1atN608IIYQQ1kvQYcrb25u8efNiZ2dH3rx58fb2jnWfqVKlwt/f//FrwIABJlRqkDAlhBBCxF5kpNUVPCvBhilvb2+8vLw4d+4cWmvOnTuHl5eXKYHqeUFBQRQpUoTjx48D0Lp1a2bNmgVA165dcXFxoUSJEgwdOvTx9/j6+lK5cmWcnZ2pUKECQUFBDBkyhEWLFlGmTBkWLVpkep1CCCFEYrd/P5QsaXUVz7KppRGe9sUXX+Dv7//S7bt27eLBgwfPfHbv3j06der0OOg8r0yZMkyePPmVx71//z5lypR53B44cCAtW7Zk6tSptG/fns8//5xbt27RpUsXAEaPHk2mTJmIiIigVq1aHDx4kKJFi9KyZUsWLVqEq6srwcHBpE6dmhEjRuDn58fUqVNj+ssghBBCiKdkygSOjpA5M9y48eL2bNnivyabDVOv83yQet3nMRV1me95derUYcmSJXTv3p0DBw48/nzx4sXMnDmT8PBwLl26xJEjR1BKkSNHDlxdXQFwdHSMVU1CCCFEUhYRAStXQrNm4OQEO3cayyDYCpsNU68bQcqbNy/nzr34AGcnJyc2bdpkej2RkZEcPXqUVKlScfPmTXLlysWZM2eYMGECvr6+ZMyYkfbt2xMaGorWWm7/F0IIIUzy66/QsSNs3Ag1a9pWkIIYzJlSSs1WSl1VSh16yfYmSqmDSil/pZSfUqqq+WW+aPTo0aROnfqZz1KnTs3o0aPj5HiTJk2iWLFiLFiwgI4dOxIWFkZwcDBp0qQhffr0XLlyhb/++guAokWLcvHiRXx9fQEICQkhPDycdOnSERISEif1CSGEEImN1sbXdu3gjz+MIGWLYjIBfS5Q/xXb1wPOWusyQEfgJxPqei1PT09mzpyJk5MTSimcnJyYOXMmnp6eseo3as5U1GvAgAGcOHGCn376ie+++45q1apRvXp1Ro0ahbOzM2XLlqVEiRJ07NiRKlWqAJA8eXIWLVpEz549cXZ2pk6dOoSGhuLu7s6RI0dkAroQQgjxGuvXQ6VKcOsW2NvD++9bXdHLKR0V+161k1J5gT+11q+cP6+UqgTM1loXe12fLi4u2s/P75nPjh49SrFir/1WkQjIz1oIIcSrbNsGX3wBy5dD7txWVwNKqb1aa5fotpmyNIJSqplS6hiwCmN06mX7eT26FOh37do1Mw4thBBCiETiwQPYsMF4X7Uq7NljG0HqdUwJU1rr5VrrokBTYOQr9puptXbRWrtkzZrVjEMLIYQQIpEYNgzq14eo+8vsEshqmKbezae13qKUKqCUyqK1vm5m30IIIYRInLQ27tAbOBCqVDGWP0hIYp35lFIF1aN1AJRS5YDkQDTLaAkhhBBCPGvRIvjgAwgPNxbjtOWJ5i/z2pEppdQCoCaQRSkVCAwFHAC01j8CHwDtlFJhwH2gpY7JrHYhhBBCJHkhIXDtGty5AxkyWF3N23ltmNJat37N9vHAeNMqEkIIIUSiFhwMAQFQtix07gwdOhjLHyRUCWRqV/xavnw5SimOHTsW7fb27duzdOnSeKnl7NmzpEqVijJlyuDs7EzlypUfP3D5TaVNmzbaz69du4abmxtly5Zl69atMe5v8uTJ3Lt377X9CyGEEE/r2BEaNoSof0LeNkh5e3uTN29e7OzsyJs3L97e3uYV+QYSbJjKnt2YrPb8K3v22Pe9YMECqlatysKFC2PfmQkKFCiAv78/Bw4c4JNPPmHMmDGm9r9+/XqKFi3K/v37qVatWoy+JyIi4oUwJYQQQsTEuHEwfz489yCTN+Lt7Y2Xlxfnzp1Da825c+fw8vKyJFAl2DB15cqbfR5Td+7cYfv27fz888+Pw5TWmh49elC8eHEaNWrE1atXH+8/YsQIXF1dKVmyJF5eXkRNF6tZsya9e/emevXqFCtWDF9fX5o3b06hQoUYPHgwAHfv3qVRo0Y4OztTsmTJGK2KHhwcTMaMGQFj1KpatWqUK1eOcuXKsWPHDgAuXbpE9erVKVOmDCVLlnxmtGnQoEE4OztTsWJFrly5gr+/P/369WP16tWUKVOG+/fvs2DBAkqVKkXJkiXp37//4+9NmzYtQ4YMwc3NjdGjR3Px4kXc3d1xd3d/af9CCCEEwNSpMGCA8b5gQahVK3b9DRo06IX/0N+7d49BgwbFruO3obW25FW+fHn9vCNHjjzTrlFD6zlzjPcPHxrtefOMtnEjZfQvrbW+ds3Y38fHaF+69MLhojVv3jzdsWNHrbXWlSpV0nv37tXLli3TtWvX1uHh4fq///7T6dOn10uWLNFaa33jxo3H39umTRvt8+iANWrU0P369dNaaz158mSdI0cOffHiRR0aGqpz5sypr1+/rpcuXao7d+78+Ptv3779Qj1nzpzRKVOm1M7Ozjp//vw6e/bs+ty5c1prre/evavv37+vtdb6xIkTOurXdMKECXrUqFFaa63Dw8N1cHDwo18zHtfXt29fPXLkSK211nPmzNHdu3fXWmv933//6dy5c+urV6/qsLAw7e7urpcvX/74+xctWvS4NicnJ33t2rXH7Zf1H53nf9ZCCCESt549tW7cWOuwMHP6U0pp4IWXUsqcAzwH8NMvyTQJdmQqrixYsIBWrVoB0KpVKxYsWMCWLVto3bo19vb2vPvuu7z33nuP99+4cSNubm6UKlWKDRs2cPjw4cfbPDw8AChVqhQlSpQgR44cpEiRgvz583PhwgVKlSrFunXr6N+/P1u3biV9+vTR1hR1mS8gIIDJkyfj5eUFQFhYGF26dKFUqVJ8+OGHHDlyBABXV1fmzJnDsGHD+Pfff0mXLh1gPDPw/Uf3nJYvX56zZ8++cCxfX19q1qxJ1qxZSZYsGZ6enmzZsgUAe3t7Pvjgg5f+2sWkfyGEEEnH1atPFuCcOBFWrIBkJqxwGRQURIoUKaLdlidPntgf4A2Zumin2TZtevLeweHZ9utkyfLs/jGZS3Xjxg02bNjAoUOHUEoRERGBUopmzZrxaCmtZ4SGhtKtWzf8/PzInTs3w4YNIzQ09PH2qB+0nZ3dMz90Ozs7wsPDKVy4MHv37mX16tUMHDiQunXrUq9ePT799FPAuIRYunTpZ47p4eFBhw4dAJg0aRLZsmXjwIEDREZGkjJlSgCqV6/Oli1bWLVqFW3btqVv3760a9cOBweHx+dhb29PeHj4C+ekX7GqRcqUKbF/xSzBmPQvhBAiaYiMhHr1IHly2LXLnBAFcPLkSTw8PHjw4AHJkyfn4cOHj7elTp2a0aNHm3OgNyAjU09ZunQp7dq149y5c5w9e5YLFy6QL18+MmXKxMKFC4mIiODSpUts3LgR4HFwypIlC3fu3HnjO/wuXrxI6tSpadOmDV9++SX79u3Dzc0Nf39//P39H49sPW3btm0UKFAAMJJ5jhw5sLOzY968eURERABw7tw53nnnHbp06UKnTp3Yt29fjGtyc3Nj8+bNXL9+nYiICBYsWECNGjWi3TddunSEhIS80TkLIYRIGuzsjNGo6dONG8TMsH79etzc3Lh27RobN25k9uzZODk5oZTCycmJmTNn4unpac7B3oBNj0y9SrZs0U82z5bt7ftcsGABA6Jmxz3ywQcfcPToUQoVKkSpUqUoXLjw43CRIUOGx5fZ8ubNi6ur6xsd799//6Vv377Y2dnh4ODA//73v2j3CwgIoEyZMmitSZ48OT/99BMA3bp144MPPmDJkiW4u7uTJk0aADZt2sS3336Lg4MDadOm5ddff41xTTly5GDs2LG4u7ujtaZhw4Y0adIk2n29vLxo0KABOXLkeBwwhRBCJF2RkTB8OBQoAO3awVP3J8WK1pqpU6fSu3dvihUrho+PD/ny5QOwJDw9T73qsk5ccnFx0X5+fs98dvToUYoVK2ZJPSJ+yc9aCCESn/Bw49Je4cLwkvGBN/bw4UN69OjBrFmz8PDwYP78+Y/nAscnpdRerbVLdNsS7MiUEEIIIWzD2bOQOTOkSwd//gmPpvDG2rVr12jRogVbtmzhq6++YuTIkdjZ2d4MJQlTQgghhHhrISFQsSLUrw9z50KqVOb0e/DgQZo0acLly5f57bffaN36lU+3s5TNhSmtdbR3zonEw6pLy0IIIcyXLh18+60RqMyyYsUK2rRpQ/r06dmyZcsbz0mObzY1VpYyZUpu3Lgh/9gmYlprbty48XgZByGEEAlPWBj06mUseQDQti0UKhT7frXWjB49mmbNmlGiRAl8fX1tPkiBjY1M5cqVi8DAQK5du2Z1KSIOpUyZkly5clldhhBCiLcUEgKrV8O775o3InXv3j06derEwoUL8fT0ZNasWaQy65phHLOpMOXg4PD4VkchhBBC2JbTpyFvXsiUCfbvNy7xmSEwMJCmTZuyb98+xo0bR79+/RLUlB+buswnhBBCCNsUEAClSxvzo8C8ILVr1y5cXV05fvw4K1eupH///gkqSIGEKSGEEELEQP78MHiwMT/KLPPmzaNmzZqkTp2aXbt20bhxY/M6j0cSpoQQQggRrbt3oWtX+O8/45EwAwYY86RiKyIign79+tGuXTsqV67Mnj17KFGiROw7toiEKSGEEEJE6/x5WLAANm0yr8/g4GCaNGnCt99+S7du3Vi7di2ZM2c27wAWsKkJ6EIIIYSw3oULkDs3FCtmzJUyK+ucOnUKDw8PTp48yfTp0+natas5HVtMRqaEEEII8djmzVCwIPj4GG2zgtSGDRuoUKECV65c4e+//040QQokTAkhhBDiKRUrwuefQ7Vq5vSntWbatGnUrVuXd999F19fX9zd3c3p3EZImBJCCCGSuJs3jQB1/z6kSAHffAMZM8a+34cPH9K1a1d69OhBw4YN2bFjB/nz5499xzbmtWFKKTVbKXVVKXXoJds9lVIHH712KKWczS9TCCGEEHFl926YORP27DGvz+vXr1O3bl1mzJjBgAEDWL58OY6OjuYdwIbEZAL6XGAq8OtLtp8BamitbymlGgAzATdzyhNCCCFEXLl6Fd55Bxo0MFY3z5HDnH4PHTqEh4cHFy9eZP78+Xh6eprTsY167ciU1noLcPMV23dorW89au4C5KFrQgghhI2bPx8KFIBDj647mRWkVq5cSaVKlQgNDWXLli2JPkiB+XOmOgF/vWyjUspLKeWnlPKThxkLIYQQ1qldGzp2NAKVGbTWjBkzhmbNmlG0aFF8fX2pUKGCOZ3bONPClFLKHSNM9X/ZPlrrmVprF621S9asWc06tBBCCCFiIDAQRowArSF7dvj+e0iVKvb93r9/H09PTwYNGkTr1q3ZsmULOXPmjH3HCYQpYUopVRr4CWiitb5hRp9CCCGEMNeSJTBhgrEQp1n+++8/qlevzsKFCxk7dizz588nlRkJLQGJ9QroSqk8wO9AW631idiXJIQQQggz3b4NGTLAF1/ABx9Anjzm9Ltnzx6aNm1KSEgIK1euTLAPKo6tmCyNsADYCRRRSgUqpToppT5TSn32aJchQGZgulLKXynlF4f1CiGEEOINjB0Lzs5w/brxsGKzgtT8+fOpXr06KVOmZOfOnUk2SEEMRqa01q1fs70z0Nm0ioQQQghhmrp1jSCVIYM5/UVERDBo0CDGjx9PzZo1WbJkCVmyZDGn8wRKHnQshBBCJDLHjsH27dCpE5Qvb7zMEBwcjKenJ3/++SefffYZU6ZMwcHBwZzOEzB5nIwQQgiRyEyYAIMHQ3CweX0GBARQqVIl/vrrL6ZNm8b//vc/CVKPSJgSQgghEgGt4e5d4/2UKcajYcx6esvGjRupUKECly9f5u+//6Zbt27mdJxISJgSQgghEoGuXaF+fXj4EFKnhty5zen3f//7H3Xr1iV79uzs2bOH9957z5yOExGZMyWEEEIkAjVrQq5cYNaVt7CwMHr16sWPP/7I+++/j7e3d6J9UHFsSZgSQgghEqg9e+DWLahXD1q1Mq/f69ev8+GHH7Jp0yb69evHmDFjsLe3N+8AiYyEKSGEECIB0hp69zbmSdWpA3YmTdw5fPgwjRs35uLFi8ybN482bdqY03EiJmFKCCGESEAiIiAy0rict3gxpEhhXpD6448/+Pjjj0mXLh2bN2/Gzc3NnI4TOZmALoQQQiQQERHg4QHduhkjUzlzghnrZWqtGTduHE2aNKFo0aL4+vpKkHoDMjIlhBBCJBD29uDiAjlyGI+GMcP9+/fp3Lkzv/32G61atWL27NlJ7kHFsSVhSgghhLBx//xj3KlXrBgMH25evxcvXqRp06b4+voyevRoBg4ciDIrpSUhEqaEEEIIG3b/PrRvD25u8Pvv5vXr6+tL06ZNCQ4OZsWKFTRp0sS8zpMYCVNCCCGEDQoLg2TJIFUqWLsW8uY1r+/ffvuNjh07kiNHDnbs2EGpUqXM6zwJkgnoQgghhI25fRtq1IBp04x2yZKQNm3s+42MjGTgwIF4enri5ubGnj17JEiZQEamhBBCCBvj6GiMRGXPbl6fISEheHp68scff/Dpp58yZcoUkidPbt4BkjAJU0IIIYSNWLECqleHTJngt9/M6/f06dN4eHhw7Ngxpk6dSrdu3WSiuYnkMp8QQghhAy5cgJYtYfx4c/vdtGkTFSpU4OLFi6xdu5bu3btLkDKZhCkhhBDCQhERxtfcuWHdOhg50ry+f/zxR+rUqcM777zDnj17qFWrlnmdi8ckTAkhhBAWOX8eypeHDRuMdrVqYMY0prCwMLp3707Xrl2pW7cuO3fupGDBgrHvWERLwpQQQghhkQwZjMnmZl51u3HjBvXq1WP69On07dsXHx8f0qdPb94BxAtkAroQQggRj7SG5cuNZ+w5OsLmzeaFqcOHD+Ph4UFgYCC//PIL7dq1M6dj8UoyMiWEEELEo+3b4YMPYO5co21WkPrzzz+pVKkS9+7dY/PmzRKk4pGEKSGEECIeaG18rVoV/vwTOnY0q1/N+PHj8fDwoHDhwvj6+lKxYkVzOhcx8towpZSarZS6qpQ69JLtRZVSO5VSD5RSX5pfohBCCJGw/fuvMdE8IMBoN2oEdiYMZ4SGhtKuXTsGDBjARx99xJYtW8iVK1fsOxZvJCY/yrlA/Vdsvwn0AiaYUZAQQgiR2KRObSyBEBJiXp8XL16kRo0azJ8/n1GjRrFgwQJSp05t3gFEjL12ArrWeotSKu8rtl8FriqlGplYlxBCCJGgRUYaDyhu0AAKFID9+80ZjQLw8/OjSZMmBAUFsXz5cpo2bWpOx+KtxOucKaWUl1LKTynld+3atfg8tBBCCBGv5s6Fhg1hyxajbVaQWrBgAdWqVcPBwYEdO3ZIkLIB8RqmtNYztdYuWmuXrFmzxuehhRBCiHgRNdG8XTtYvNhYiNMMkZGRDBo0iI8//hhXV1d8fX0pXbq0OZ2LWJG7+YQQQgiTbNkC7u7G3KhkyeDDD81Z+iAkJIRmzZoxZswYunTpwrp165BBCdshi3YKIYQQJgkLgxs34OZNSJfOnD7PnDmDh4cHR48eZcqUKfTo0UMeVGxjXhumlFILgJpAFqVUIDAUcADQWv+olMoO+AGOQKRS6guguNY6OM6qFkIIIWxEWBj4+kLlylCrFvj7g729OX1v3ryZDz74gIiICNasWUPt2rXN6ViYKiZ387V+zfbLgCxqIYQQIkkaMgQmToSTJyFPHvOC1IwZM+jRowcFCxbEx8eHQoUKmdOxMJ1c5hNCCCFioW9fY0HOPHnM6S8sLIzevXszbdo0GjRowIIFC+RBxTZOJqALIYQQb2jFCmjTxlhLKlMmaNHCnH5v3rxJ/fr1mTZtGl9++SV//PGHBKkEQEamhBBCiDd04YJxWS84GDJkMKfPI0eO4OHhwYULF5g7dy6ffPKJOR2LOCcjU0IIIUQM3L0Lhw8b73v0gK1bzQtSq1atomLFity5c4dNmzZJkEpgJEwJIYQQMdC+PdStC/fuGWtHJU8e+z611nz77bc0btyYQoUK4evrS6VKlWLfsYhXcplPCCGEiIHhw+HcOeOhxWYIDQ3Fy8uLefPm8dFHHzFnzhx5UHECJSNTQgghxEv89JMRogCKFzceWmyGS5cuUbNmTebNm8eIESNYuHChBKkETEamhBBCiJfYvduYbB4ebjwexgx+fn40bdqUW7dusWzZMpo3b25Ox8IyEqaEEEKIp9y4AaGhkDMnTJsGdnbmBalFixbRvn17smXLxo4dO3B2djanY2EpucwnhBBCPBIZCXXqwEcfgdbGJHMzglRkZCSDBw+mVatWuLi4sGfPHglSiYiMTAkhhBCP2NnBmDHGkgdmPUv4zp07tG3blhUrVtCpUyemT59OcjNuBRQ2Q8KUEEKIJE1rGDcO8ueHli2hfn3z+j5z5gxNmjTh8OHDfP/99/Ts2RNlVkoTNkMu8wkhhEjSwsNh9Wr45x9z+92yZQsVKlTgwoULrFmzhl69ekmQSqRkZEoIIUSS9N9/kDGjsW7UX39BmjTm9T1r1iy6detGgQIF8PHxoXDhwuZ1LmyOjEwJIYRIckJCwNUVevUy2mnTmjNHKjw8nJ49e+Ll5UXt2rXZtWuXBKkkQEamhBBCJDnp0hmLcVapYl6fN2/e5KOPPmL9+vX83//9H9988w329vbmHUDYLBmZEkIIkSSEh0P//rBvn9Hu0sVY1dwMR48exc3Nja1btzJnzhy+++47CVJJiIxMCSGESBJu34YFC4xRqXLlzOt39erVtG7dmpQpU7Jx40YqV65sXuciQZAwJYQQIlG7cAFy5YIsWcDfHzJlMqdfrTXfffcd/fr1w9nZmZUrV5InTx5zOhcJilzmE0IIkWgFBEDJkjBpktE2K0iFhobSvn17+vbtS4sWLdi2bZsEqSRMwpQQQohEK39++OILaNHCvD4vX76Mu7s7v/76K8OHD2fRokWkMXNdBZHgSJgSQgiRqISGwv/9H1y5Yix3MHw4mDVotG/fPlxdXTl48CBLly5lyJAhshCnkDlTQgghEq7s2Y3QFJ0yZaBdO/OOtWjRIjp06EDWrFnZvn07ZcqUMa9zkaC9dmRKKTVbKXVVKXXoJduVUmqKUuqUUuqgUsrEeySEEEKIl3tZkALzglRkZCRDhgyhVatWlCtXDl9fXwlS4hkxucw3F3jVYx8bAIUevbyA/8W+rNjJnt0Y2n3+lT271ZUJIYRISO7cuUOLFi0YOXIkHTt2ZP369bzzzjtWlyVszGvDlNZ6C3DzFbs0AX7Vhl1ABqVUDrMKfBsv+5/Kq/4HI4QQImEJD4/b/s+ePUuVKlVYuXIlkyZN4qeffiJFihRxe1CRIJkxZyoncOGpduCjzy6Z0LcQQggRraZN467vrVu30rx5c8LCwli9ejX16tWLu4OJBM+Mu/miu41BR7ujUl5KKT+llN+1a9dMOLQQQoikIjISfHyejEh17x43x/npp5+oVasWmTJlYvfu3RKkxGuZEaYCgdxPtXMBF6PbUWs9U2vtorV2yZo1qwmHfnPBwQ8sOa4QQojYWbcOmjSBZcuMdoMGkC1b9Pu+7PNXCQ8P5/PPP6dLly64u7uza9cuihQp8vYFiyTDjDDlA7R7dFdfRSBIa22zl/gqVmzDiRMnrC5DCCHEa2gNq1bB0qVGu04dY2Tq6QU4L1829nv+dfnymx3r1q1bNGzYkClTptC7d29WrVpFxowZzTsZkajFZGmEBcBOoIhSKlAp1Ukp9ZlS6rNHu6wGTgOngFlAtzirNoZe9j8SR8cHXLmygXLlyjFihA+RkdFejRRCCGEjvv3WeBSM1sZd2Y0bg729ucc4duwYbm5ubNq0iZ9//pmJEyeSLJkswyhiTmltTaBwcXHRfn5+8X7cwMBAPDxGsX//dNzcZvP33x/h6OgY73UIIYR40f798PXX4O0N6dPDxYuQNSs4OMTN8dasWUOrVq1IkSIFv//+O1WqVImbA4kETym1V2vtEt22JPc4mVy5crFr1zQaNvybPXt6U65cOXbv9rW6LCGESNIiIoyvWoO/P0TNxnj33bgJUlprJk6cSKNGjcibNy++vr4SpMRbS3JhCiB5cntWrarP1q1rePgwkkqV7tKs2XoiIyOtLk0IIZKU8HBo2BAGDDDa5crBmTPg6hp3x3zw4AEdO3akT58+NGvWjO3bt5PHrIf3iSQpSYapKFWqVGHbtr28+24KVqz4iYYNG3JFVvYUQog4999/xtdkyaBIkWcfRBxXl/QALl++jLu7O3PnzmXYsGEsXryYNGnSxN0BRZKQpMMUQJ48GTl/viI//liTzZs3U6TIl4wZs9fqsoQQItGaORPy54ezZ432pEnQs2fcH3f//v24urri7+/PkiVLGDp0KHZ2Sf6fQWEC+V0E2NkpPv30U3x9/QgL+5xBg8Lo27c/Dx8+tLo0IYRIFI4cgdOnjfcNG8JXX4WArUcAACAASURBVEGmTPF3/CVLllClShWUUmzfvp0WT6+vIEQsSZh6SsmSJbhwoQTt2q1kwoRvqFy5Fhs3nrO6LCGESNDu3oVKlWDYMKOdKxcMHQrxcSN1ZGQkQ4cO5aOPPqJs2bL4+vpStmzZuD+wSFIkTD0nU6ZU/PLLWJYtW8ahQ614770MTJ++1OqyhBAiQTl9GiZPNt6nSQOLF8PEifFbw927d/nwww8ZMWIEHTp0YMOGDWR7m6XRhXgNCVMv0bx5czZubEqBAjPp3v1D2rdvT0jIHavLEkKIBGHxYhg4EAIDjXa9epAlS/wd/9y5c1SpUoUVK1YwceJEfv75Z1KkSBF/BYgkRcLUK1SqlJNjx3ozZMgQfv11L1myXGDp0iNWlyWEEDbn9m3jwcMbNhjtnj0hIMC4pBfftm3bhqurK2fPnmXVqlX07t0bpVT8FyKSDAlTr5EsWTKGDx/O1Km/oPVDWreux+TJk7Fq5XghhLAlUcvzpUoFa9YYC26CcWnv3XfjpwZvb2/y5s2LnZ0dmTNnpkaNGmTIkIHdu3dTv379+ClCJGkSpmKoW7dyXLyYi0aNytO7d2+KF5/NkSPXrS5LCCEsM3kyVKliBKoUKYw79v7v/+K3Bm9vb7y8vDh37hxaa27evAnAl19+SZEiReK3GJFkSZh6A1myZGb58uUMGuTNsWOeVKw4hfXr11tdlhBCxJubN41Vy8F4qHzBgsbdemAEqvg2aNAg7t2798xnkZGRjBkzJv6LEUmWhKk3pJRi1KiPWbnyLDlzLqNOnTp07fot9+6FWV2aEELEqZMnIV8+mDfPaLdubbxPl86aeg4dOsS5c9EvX3P+/Pl4rkYkZRKm3pKHR1H8/PbQvn1XfvyxBXnybOTMmTNWlyWEEKYKCQHfR8+CL1gQPvsMKlSwtqZz587Rvn17Spcu/dKJ5fKsPRGfJEzFQpo0aZg9expffHGJ0NDxlClThoULF1ldlhBCmKZdO/DwgIcPQSkYPx5KlLCmluvXr9O7d28KFy7MwoUL+fLLL/nxxx9JnTr1M/ulTp2a0aNHW1OkSJIkTJlg0qTKHDr0M8WLF6d1a3+KFNnK7dt3rS5LCCHe2P37MGUKBAUZ7aFDwccHkie3rqY7d+4wcuRI8ufPz5QpU2jbti2nTp3im2++wcvLi5kzZ+Lk5IRSCicnJ2bOnImnp6d1BYskR1l1i7+Li4v28/Oz5NhxJSwsjPfe28q2bZcoUmQUixYtxNnZ2eqyhBAixvbvh3LlYPZs6NDB2loePnzIrFmzGDlyJFeuXKFZs2aMHj2aYsWKWVuYSJKUUnu11i7RbZORKRM5ODiwdet7/PNPDoKDg6hQoQmtWm0kMlLWpBJC2K6ff4ZvvzXely0L//5rbZCKjIzkt99+o1ixYvTo0YOiRYuyc+dOfv/9dwlSwiZJmIoDtWu/x4EDB3ByGsOiRW7Uq9eFGzduWF2WEEI89vRFic2bYe3aJ5+VLGlVTZo1a9ZQvnx5PD09SZcuHatXr2bjxo1UrFjRmqKEiAEJU3Eka9asHDvWmv79l7FlyzycnZ357bddVpclhBDs2mUEprNnjfaPP8I//xgTzK2ye/du3nvvPRo0aEBQUBDe3t7s27ePBg0ayKNghM2TMBWH7OwU48a1ZdeuXUBtPD1daN16PuFRK94JIUQ8iYw0np8HkDs3ODpC1IB56tTWBaljx47RvHlzKlasyJEjR/jhhx84duwYH3/8MXZ28k+USBjkd2o8KFu2LH5+U3F2XsXChZ2pWbMmZ89Gv9CcEEKYTWuoWhU+/dRo58wJO3dC+fLW1RQYGEjnzp0pUaIE69atY/jw4QQEBNCjRw+SW3nroBBvQcJUPMmePS3+/k3w9p7NgQNHKVToNP3777C6LCFEIqW1EZjAGHXy9ITmza2tCeDmzZv069ePQoUKMW/ePHr16kVAQABDhgwhbdq0VpcnxFuRpREssGvXWWrXvsPdu1/z6afZmDhx4guLzgkhRGz8/DN07gy7d1u/YjnAvXv3mDJlCuPHjycoKIi2bdsyfPhw8ubNa3VpQsRIrJdGUErVV0odV0qdUkoNiGa7k1JqvVLqoFJqk1IqV2yLTswqVszLtWuF6devMDNmzKBYsb4sX37S6rKEEAmY1rB+vRGeAFq2hDlzjKUOrBQWFsbMmTMpVKgQAwcOpGrVqhw4cIBffvlFgpRINF47MqWUsgdOAHWAQMAXaK21PvLUPkuAP7XWvyil3gM6aK3bvqrfpDwy9bS//vqHxo3zo/UZpk07xaeffip3rggh3tjDh1CgALi6wu+/W12NsczB0qVLGTx4MCdOnKBy5cqMHz+eqlWrWl2aEG8ltiNTFYBTWuvTWuuHwEKgyXP7FAfWP3q/MZrt4iUaNKiDv386qlWbQ9euXWnatA0BAbesLksIkQD4+UHXrsadesmTw19/wW+/WV0VrF+/ngoVKvDRRx/h4OCAj48P27ZtkyAlEq2YhKmcwIWn2oGPPnvaAeCDR++bAemUUpmf70gp5aWU8lNK+V27du1t6k2USpZ8hw0b5jFhwgT++KM2RYqE8Pff260uSwhho6IuKBw/boxCnT5ttEuWhJQpratr79691K1bl9q1a3P16lXmzp3LgQMHaNy4sYy4i0QtJmEquj8Bz18b/BKooZTaD9QA/gNeWExJaz1Ta+2itXbJmjXrGxebmNnZ2dGnTx/mznUjY8bZNGhQnREjRhAREWF1aUIIGxEUBB4eMHeu0W7VCgICoGBBS8vi5MmTtGzZEhcXF/bt28fEiRM5fvw4n3zyCfb29tYWJ0Q8iEmYCgRyP9XOBVx8eget9UWtdXOtdVlg0KPPgkyrMglp1644p0/34eOPP2bo0OVkybIfP7+Lr/9GIUSiFRJifHV0hNBQY34UgL09WLmawKVLl+jatSvFixdn1apVfP3115w+fZrevXuT0sohMiHiWUzClC9QSCmVTymVHGgF+Dy9g1Iqi1Iqqq+BwGxzy0xa0qVLx7x58/j880kEBTlSp04NVq5caXVZQggLjB0LRYvCvXvGelFr1z5ZfNMqQUFBDBo0iIIFC/LTTz/x6aefEhAQwIgRI3B0dLS2OCEs8NowpbUOB3oAa4GjwGKt9WGl1AillMej3WoCx5VSJ4BswOg4qjdJmTy5JocPK/Lnd6Rp06bUrLmQ27dDrS5LCBHHTp2C4GDjffXq8MknEHXF38qpR6GhoXz33Xfkz5+fMWPG0KRJE44dO8bUqVPJli2bdYUJYTFZtDMBePjwIW3bzmTx4h7kyjWIv/9uQ7FixawuSwgRBwIDIV8++PprGDLE6moM4eHh/PrrrwwdOpTAwEDq1avH2LFjKWv1IlZCxKNYL9oprJU8eXIWLerBpEk7CA2dRfny5ZkwwZvISGuCsBDCXBcuPFkbKlcu+PFH8PKytiYw1opasWIFpUuXplOnTrz77rts2LCBNWvWSJAS4ikSphKQL76ozMGDB6hQoT59+9akeHEfbkc9Bl4IkWANGQIdO8Ldu0a7UyfInt3amrZs2UKVKlVo1qwZkZGRLFu2jF27duHu7m5tYULYIAlTCUyOHDlYt24pDRqc4NSpUZQpU4adUU8zFUIkCDduQJ8+cPas0R45Eg4cgDRpLC0LgAMHDtCwYUNq1KjB+fPnmTVrFocOHaJ58+ayVpQQLyFhKgFKlsyO1avd2b59KnZ2dlSpspY6dTYRHi5rUgmREISGwowZsGGD0c6VC5ycrK3pzJkztGnThrJly7Jr1y6++eYbTp48SefOnUmWLJm1xQlh4yRMJWBubm7s27efPHlqs27daerWrcPFi7ImlRC2aNw46NzZeJ8zJ5w/b1zas9rVq1fp1asXRYoUYdmyZfTv35+AgAD69u1LqlSprC5PiARBwlQClyFDek6frsLMmXbs3r2bEiUaMHToHqvLEkJgrA0V5e5duHPnyRIHmTJZU1OU4OBghg4dSv78+Zk+fTodOnTg1KlTjB07lowZM1pbnBAJjISpRMDOTtGlS3v27t2Lnd0QRowowmeffcWDBw+sLk2IJGvnTsid2/gKMGIELFxorFpupQcPHvD9999ToEABRowYQYMGDTh8+DAzZswgZ87nH7sqhIgJCVOJSNGiRQkIaETLlj8xY8ZYKlasyI4dp6wuS4gk4949OHPGeF+qFNSvbzwCBqxdbBMgIiKCefPmUbRoUb744gtKly7Nnj17WLJkCUWKFLG2OCESOAlTiUyGDClZuLAPPj4+BASUpEqV7AwZ8idWLc4qRFLi7g4ffwxaG8/M8/aGEiWsrUlrzapVqyhbtizt2rUjU6ZM/P3336xbtw5XV1drixMikZAwlUg1btyYDRu+IWfOvxk5sjmenp4ERz2fQghhigcPYP58iIw02kOGwDffWD8KFWXHjh3UqFGD999/n/v377Nw4UJ8fX2pU6eOLHMghIkkTCViLi45OHeuCaNGDWXRIh9y5NjNnDmHrS5LiETDxwfatoV164x2o0ZQrZq1NQEcPnyYJk2aUKVKFU6ePMn06dM5cuQILVu2xM5O/toXwmzypyqRs7e3Z9CgQXh7b+HBg6J07tyPb775hsio/0oLIWIsMhJ++eXJo1+aN4eNG6FOHWvrinL+/Hk6dOhA6dKl2bRpE6NGjeLUqVN07doVBwcHq8sTItGSMJVEtGpVjv/+S0fz5qnp378/5ct/zcGDV60uS4gERSmYOtWYCwXGnXk1a1p/We/69ev06dOHwoULs2DBAnr37s3p06cZNGgQaWxhWXUhEjkJU0lItmwZWLx4MVOmzMbfvzcVKmxnzZo1VpclhE3btAlq1TLu1FMKVq2CpUutrspw9+5dRo0aRYECBZg8eTIff/wxJ06cYMKECWTOnNnq8oRIMiRMJTFKKXr27MDKlbfJl+8HGjRoQK9eg7hz56HVpQlhM7SGh4/+SNjbw6VLxorlAO+8Y/1IVFhYGNOnT6dAgQJ8/fXXuLu7c/DgQWbPnk2ePHmsLU6IJEjCVBLl4VGQfftW0a1bN374oTTZsx/j2DFZk0qI+/ehQgUYO9ZoV6sGhw5B0aLW1gUQGRnJwoULKVasGN27d6dw4cJs376dFStWUMLqNRiESMIkTCVhqVKlYtq0aQwYkJfISG9cXcsyf/58q8sSIt5pDSdPGu9TpYKqVZ8NT1bfAKe15u+//8bFxYXWrVuTJk0aVq1axebNm6lcubK1xQkhJEwJGDvWjePHe1C2bFnatp1GwYJbuHgxxOqyhIg3o0dD6dLG5TyASZOgZUtra4qyZ88eatWqRb169bh16xbz5s1j//79NGzYUNaKEsJGSJgSAOTOnZsNGzbQsOHXBATkoGrV6uzdu9fqsoSIM9u2wblzxvvWreG776x/+PDTjh8/TosWLXBzc+PQoUN8//33HDt2jDZt2shaUULYGPkTKR5LliwZq1Y15J9/rhIWdp2KFSvTtu0qwsNlTSqRuFy/DrVrGwEKoEAB6NYNUqSwti6A//77Dy8vL0qUKMHatWsZOnQoAQEB9OrVixS2UKAQ4gUSpsQLateuwoEDByhbdiTz5zfCzW0YV6/KmlQiYdu3DyZMMN5nyQKrV8O4cdbW9LRbt27Rv39/ChYsyNy5c+nevTsBAQEMGzaMdOnSWV2eEOIVJEyJaGXKlIldu/rSvfufHDr0Dc7OzqxcudHqsoR4a7//boSnW7eM9nvvQerU1tYEcO/ePcaPH0/+/Pn59ttvadGiBcePH+f777/nnXfesbo8IUQMSJgSL2Vnp5g69X38/HxxdCxM06bFqFt3BWFhYVaXJsQzsmc31n56/pUyJezcaezTvz8EBEDGjNbWGiU8PJxZs2ZRqFAhBgwYQOXKlfH392fevHnky5fP6vKEEG8gRmFKKVVfKXVcKXVKKTUgmu15lFIblVL7lVIHlVINzS9VWKVUqVJs3/4XxYsf5Z9/vqJq1aqcPn3a6rKEeOzKleg/f/AATj1aPi1dOkifPv5qehmtNUuXLqVEiRJ4eXnh5OTE5s2bWbVqFaVLl7a6PCHEW3htmFJK2QPTgAZAcaC1Uqr4c7sNBhZrrcsCrYDpZhcqrJUlS2oOH3ZnyZIRnDhxgmLFFtGz5w6ryxKCe/devb1t2/ipIyY2bNiAm5sbH374IcmSJWPFihVs376d6tWrW12aECIWYjIyVQE4pbU+rbV+CCwEmjy3jwYcH71PD1w0r0RhS1q0aMGePf6kSNGQqVP307FjR+7evWt1WSKJ0Nq4VBdl5EjjEp+t279/P/Xq1aNWrVpcvnyZ2bNnc/DgQZo0aSJrRQmRCMQkTOUELjzVDnz02dOGAW2UUoHAaqBndB0ppbyUUn5KKb9r1669RbnCFhQq5MTVqyUYMOAac+fOpVSppixadNzqskQi9OAB7NjxZPRpxgwoWPDJc/KqV4d+/ayr73VOnTpF69atKVeuHH5+fkyYMIETJ07QoUMH7O3trS5PCGGSmISp6P7bpJ9rtwbmaq1zAQ2BeUqpF/rWWs/UWrtorV2yZs365tUKm5EyZTLGjh3G+vXruXTpS1q1cuS776ah9fO/NYSIubt3YcWKJ2Fp82aoUsUIVAB16xqBKmqlgBo1YPBga2p9lcuXL9O9e3eKFSuGj48PgwYN4vTp0/Tp04eUKVNaXZ4QwmQxCVOBQO6n2rl48TJeJ2AxgNZ6J5ASyGJGgcK2ubu74+/vQuXKk/nyyx54eHhw9uwNq8sSCURoKMyaBbt3G+1r16BZM/DxMdqVK8Py5eDiYrTz5wcvrxfvyMuWLfr+X/Z5XAkKCmLw4MEUKFCAmTNn0qVLF06dOsWoUaNIbwuz34UQcSImYcoXKKSUyqeUSo4xwdznuX3OA7UAlFLFMMKUXMdLIooUycy2beOYMmUKf/2ViQIFwpg/f5fVZQkbpLWx1tPChUbb3h6++AKWLTPaTk7GUgadOxvttGmhaVPIkOHV/V6+bPT9/Ovy5bg7l6eFhoYyceJEChQowOjRo2ncuDFHjx5l+vTp5MiRI36KEEJY5rVhSmsdDvQA1gJHMe7aO6yUGqGU8ni0Wx+gi1LqALAAaK/lek+SopSiZ8+eeHt/Rdq0O2jbthqDBw8mPDzc6tKExUaOhGHDjPdKwYIFsH690XZwgOPHYfz4J9srVjTWh0oIIiIimDNnDoULF6ZPnz6UL18ePz8/Fi5cSMGCBa0uTwgRT5RVmcfFxUX7+flZcmwRt+7evUuvXr2YPXsRWbIsYfnyUlStmsvqskQc0toIQgBjxoCvr3F5DqBdO4iIAG9vo/3ggW08Ay82tNb4+Pjw1VdfceTIEVxdXRk3bhzvvfee1aUJIeKIUmqv1tolum2yArowXZo0afj5558ZPnw5169XpH79z1iyZInVZQkT3bnz5P3kycbluchHz8NOkcIYWYr6f9qvvz4JUlHbE7KtW7dStWpVmjZtSnh4OEuWLGH37t0SpIRIwiRMiTgzZEgd/P1vU7LkdT766CMaNZrI9euvWWFR2Byt4exZY0QJ4JdfjJXEL10y2oUKwfvvG3fiAfTpY1zKS2zLJ/3777+8//77VK9enTNnzjBjxgwOHTpEixYtZK0oIZI4CVMiTjk752Pr1q307DmK1au7Urjw7/z7779WlyVeISwM9uyBq1eN9t9/Q758sOvRPQWurvD112D36G+PRo1g+vQnyxUkNmfPnqVdu3Y4Ozuzbds2xo4dy6lTp/Dy8sLBwcHq8oQQNkDClIhzDg4OTJkyiG++OYKDw1hcXV2ZPHkmkZFyj4ItePgQVq2CQ4eM9tmz4OZmrPcExvupU40RKIDixY0J5fG97EB8u3btGp9//jmFCxdmyZIl9O3bl9OnTzNgwABSp05tdXlCCBsiE9BFvLp69SqffNKeNWs+I2fOTBw4UJzMmTNZXVaSorUxj+mdd6BBA+PyXfr00L07fPedsX35cqhWDZLi2rohISFMnDiRCRMmcO/ePTp27MjQoUPJlUtuohAiKZMJ6MJmvPPOO/z55580aeLIpUuLKVPGmS1btlhdVqI3Ywb88IPxXiljrac5c4x2ihTGCuMjRjzZ3rx50gtSDx484IcffqBAgQIMGzaMunXrcvjwYWbNmiVBSgjxShKmRLyzt7djxYqa7NnzCSlTpqRmzWHUrLmJ0FBZk8os//sffPbZk/ZffxmX8qJs3Phk4UyAcuUgTZr4q8+WREZGMn/+fIoWLUqvXr0oWbIku3fvZtmyZRQtWtTq8oQQCYCEKWGZ8uXLs2/fPooX78PmzTmoVasRFy5ceP03CuDJUgQAP/8MFSo8WY7g0iU4efJJe/FiWLPmyf7Zsz+ZQJ5Uaa1ZvXo1ZcuWpW3btmTIkIE1a9awfv16KlSoYHV5QogEJIn/dSqsli5dOg4dasSPPx7g4MEdlC5djhEjtlpdlk26exeiFpRfvBiyZDGeZQfGY1dy5oTgYKM9YoSxynjUHfvJk8d/vbZs586d1KxZk0aNGnHnzh1+++039u7dS7169WSZAyHEG5MwJWzCp59+xP79+3F0/JKhQ6vRosV47t+/b3VZlgoMhKAg4/2GDcYk8T17jHbBgtCihfGgYICWLY1J4/Is3Vc7cuQITZs2pXLlyhw7doypU6dy9OhRWrdujV1SH6oTQrw1+dtD2IyCBQty+HBv3n9/PsuWDaBChQrs23fE6rLiRWQk7N8PZ84Y7ePHIXdu+P13o12qFPTv/2Q5gnLlYOZMYx/xehcuXKBjx46UKlWKDRs2MHLkSAICAujevTvJZdhOCBFLEqaETUmbNjl//NGGNWvWcOmSA+XLZ6Br179JbM/N1hrWroXt2432gwfGek4zZhjtQoWMu++qVzfaWbPC6NFQoIA19SZUN27c4Msvv6RQoUJ4e3vz+eefc/r0aQYPHkzatGmtLk8IkUjIOlPCZv3771Vq1z7O1avtadGiHDNnziRjxoxWl/XWfv8d7t8HT0+jXbAglC79ZPRpzRpjBCpnTutqTMi8vb0ZNGgQ58+fJ1euXFSsWJG1a9cSEhJCu3btGD58OE5OTlaXKYRIoF61zpSEKWHTIiMj+e677/jqq69Ik2Y848bV4bPPSlldVowsWgSHDz9Zv6lhQ7h+/cm8pyNHjMt0ifUxLPHJ29ubLl26vDDPrly5cvzyyy+ULFnSosqEEImFhCmR4G3Y4Ee9elmJiFjGiBH3GDhwIPb29laX9YwlS8Db25gIrhT07m081+7QIaN9/TpkzAg2VrblwsPDCQ4OJigoiKCgoGfex7R9+/btaPt2cnLi7Nmz8XtCQohE6VVhKll8FyPE23jvPRfOnAmmT5/9fP31fP744xjTpn2Li0uOeK1Da+NlZwd//gn9+hnznjJmNO68u3TJ+JohA3z7LUya9OR7s2SJ11LjnNaaBw8evFHwie6ze/fuvfZYDg4OpE+f/plX/vz5H7//IWp59+ecP3/e7NMWQogXSJgSCUauXI4sXPgrDRrUplOnwlSsGMyyZX40adI4zo4ZddUoVSrYtg2aNjVWE3d1NQJU/vxw+7bxvnNn4xUlmQ3/6YqMjOTOnTtvNQr0dDssLOy1x0qTJg3p06fH0dHxcfjJkyfPM+3ntz/fTpky5SuP4ePjw7lz5174PE+ePG/9aySEEDFlw3/dC/EipRTt239C9uyn6dVrDE2b/kyPHj0ZPfobHB1f/Q9uTFy5YnzNlg1OnICSJY1n2Hl6GsGpcWMjWAFUqWKMTsW38PDwtx4FimoHBwe/9g5JOzu7ZwKNo6MjOXLkoGjRoq8NP1FtR0dHksVDqhw9ejReXl7PjHKlTp2a0aNHx/mxhRBC5kyJBOvBgwcMGDCAyZNDSJmyL+vW2VGlSqEYf7/WxnwmMO6iu3cPHB1h4EAYORIiIuDrr40FMZ2dY1+v1prQ0NBYzQ2K6WWx5MmTv3BZ7E1GgtKnT0+aNGkS1GrgT9/NlydPHkaPHo1n1K2TQggRSzIBXSRqw4fvYezYCzx4UBXI9sL2bNng8mXj/bZtEBICDRoYbScnqFjRuPMOYO5ccHExRqSeFnVZLDaXxIKDg2N0WSxt2rRvHHyeb6dIkeLtf0GFEEK8QMKUSPQuXrxIzpzvvnT79es3CAoK4uOPs3DjhuK77zYSFBTE3r1psLe/SLJkF14ZhEJCQmJ0Wextw09U29HR0ebuUhRCCCF384kk4N13Xx6kALJkyQpoID9wgyZNgp7ZniJFihfCTaFChd4oDCW0y2JCCCHMIWFKJAnffz/5pWHI0dFRLosJIYR4azEKU0qp+sD3gD3wk9Z63HPbJwHuj5qpgXe01hnMLFSI2OjVq5fVJQghhEikXhumlFL2wDSgDhAI+CqlfLTWR6L20Vr3fmr/nkDZOKhVCCGEEMLm2MVgnwrAKa31aa31Q2Ah0OQV+7cGFphRnBBvItuLN/K98nMhhBDCDDG5zJcTuPBUOxBwi25HpZQTkA/Y8JLtXoAXyMrEwnxRyx8IIYQQ8SkmI1PR3Z70snvEWwFLtdYR0W3UWs/UWrtorV2yZs0a0xqFEEIIIWxWTMJUIJD7qXYu4OJL9m2FXOITQgghRBISkzDlCxRSSuVTSiXHCEw+z++klCoCZAR2mluiEEIIIYTtem2Y0lqHAz2AtcBRYLHW+rBSaoRSyuOpXVsDC7VVS6oLIYQQQlggRutMaa1XA6uf+2zIc+1h5pUlhBBCCJEwWPZsPqXUNeBcPBwqC3A9Ho5ji+Tck66kfP5J+dwhaZ+/nHvSFR/n76S1jvbuOcvCVHxRSvm97MGEiZ2ce9I8d0ja55+Uzx2S9vnLuSfNcwfrzz8mE9CFEEIIIcRLSJgSQgghhIiFpBCmZlpdgIXk3JOupHz+SfncIWmfv5x70mXp+Sf6OVNCCCGEEHEpKYxMCSGEEELEGQlTQgghhBCxsVhzlAAAB1FJREFUkGDDlFKqvlLquFLqlFJqQDTbUyilFj3avlsplfepbQMffX5cKVUvPus2QwzO/f+UUkeUUgeVUuuVUk5PbYtQSvk/er3wWKCEIAbn314pde2p8+z81LZPlFInH70+id/KYy8G5z7pqfM+oZS6/dS2BP2zV0rNVkpdVUodesl2pZSa8ujX5qBSqtxT2xL0zx1idP6ej877oFJqh1LK+altZ5VS/z762fvFX9XmiMG511RKBT31+3vIU9te+WfG1sXg3Ps+dd6HHv05z/RoW4L+uQMopXIrpTYqpY4qpQ4rpT6PZh/r/+xrrRPcC7AHAoD8/9/e3YfIddVhHP8+2lhpG2xifIl9wVqKtAnVtEU0FbUqWCPNtlRhfalUimJbK/6hIAoKFaFQBP/QUt8KVjSxRiu2tNpAKkJrqrUY0pRi043UEEFItEmMrMQ8/nHO6Ox0N3t3Z3dn7uzzgWHPnDl3OL/5zRnOnjlzL/ASYBdwUU+bm4A7a3kc+HEtX1TbnwqcV5/nxYOOaYFjvwI4rZZv7MRe7x8ddAxLEP/1wDemOXY1MFH/rqrlVYOOaSFj72l/C3DXCOX+bcAlwJMzPL4JeBAQ8GbgsVHI+xzi39iJC3hvJ/56/8/AmkHHsIixvwO4f5r6OY2ZYbzNFntP26uAHaOS9xrDWuCSWl4J/Gmaz/yBj/22rky9Cdhre8L2v4GtwFhPmzHg+7W8DXiXJNX6rbYnbe8D9tbna4tZY7f9sO1j9e5O4Owl7uNiapL7mbwH2G77kO2/A9uBKxepn4thrrF/ENiyJD1bArZ/Axw6SZMx4G4XO4EzJa2l/XkHZo/f9qM1Phixcd8g9zPp5/NiKMwx9pEa8wC2/2r7iVo+QrlG8Fk9zQY+9ts6mToL+EvX/f288MX9XxuXizU/D7y84bHDbK79v4EyY+94qaTHJe2UdPVidHCRNY3/2rrcu03SOXM8dlg17n/9avc8YEdXddtzP5uZXp+2530+ese9gYck/UHSJwbUp8X2Fkm7JD0oaV2tWza5l3QaZaLw067qkcq7ynadDcBjPQ8NfOw3utDxENI0db3neJipTZNjh1nj/kv6CHAZ8Pau6nNtH5D0OmCHpN22n12Efi6WJvHfB2yxPSnpk5QVync2PHaYzaX/48A22//pqmt77mczqmN+TiRdQZlMvbWr+vKa+1cC2yU9XVc8RsUTlOumHZW0Cfg5cAHLK/dXAY/Y7l7FGpm8SzqDMlH8jO3DvQ9Pc8iSjv22rkztB87pun82cGCmNpJOAV5GWSptcuwwa9R/Se8Gvghstj3Zqbd9oP6dAH5NmeW3yazx2z7YFfN3gEubHjvk5tL/cXqW+0cg97OZ6fVpe94bk3Qx8F1gzPbBTn1X7v8G3Eu7tjbMyvZh20dr+QFghaQ1LKPcc/Ix3+q8S1pBmUj90PbPpmky+LE/yI1l871RVtQmKF9jdDYVrutpczNTN6DfU8vrmLoBfYJ2bUBvEvsGyqbLC3rqVwGn1vIa4BnatxmzSfxru8rXADtreTWwr74Oq2p59aBjWsjYa7vXUzaeapRyX/v+WmbehPw+pm5C/d0o5H0O8Z9L2QO6saf+dGBlV/lR4MpBx7LAsb+6836nTBieq++DRmNm2G8ni70+3lksOH0E8y7gbuDrJ2kz8LHfyq/5bB+X9CngV5Rfa9xle4+kW4HHbf8C+B7wA0l7KW+y8XrsHkn3AE8Bx4GbPfWrkKHWMPbbgTOAn5Q99zxnezNwIfAtSScoq5K32X5qIIHMU8P4Py1pMyW/hyi/7sP2IUlfAX5fn+5WT10SH2oNY4eyCXWr66dJ1frcS9pC+dXWGkn7gS8DKwBs3wk8QPlVz17gGPCx+lir897RIP4vUfaF3lHH/XHblwGvAu6tdacAP7L9yyUPoA8NYn8/cKOk48C/gPH6/p92zAwghHlrEDuUfxofsv3PrkNbn/fqcuA6YLekP9a6L1D+eRiasZ/LyURERET0oa17piIiIiKGQiZTEREREX3IZCoiIiKiD5lMRURERPQhk6mIiIiIPmQyFRFDT9KZkm6q5ddI2jboPkVEdOTUCBEx9Oo1ue63vX7AXYmIeIFWnrQzIpad24Dz60n7ngEutL1e0vXA1ZQTMq4HvkY50/V1wCSwqZ6473zgm8ArKCf1+7jtp5c+jIgYRfmaLyLa4PPAs7bfCHyu57H1wIcolxH5KnDM9gbgt8BHa5tvA7fYvhT4LHDHkvQ6IpaFrExFRNs9bPsIcETS88B9tX43cHG92vxG/n95JSjX5oyIWBCZTEVE2012lU903T9B+Yx7EfCPuqoVEbHg8jVfRLTBEWDlfA60fRjYJ+kDACresJCdi4jlLZOpiBh6tg8Cj0h6Erh9Hk/xYeAGSbuAPcDYQvYvIpa3nBohIiIiog9ZmYqIiIjoQyZTEREREX3IZCoiIiKiD5lMRURERPQhk6mIiIiIPmQyFREREdGHTKYiIiIi+vBfRm+S+k45secAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "tags": [], "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "5GQLL-VobVQy" }, "source": [ "## Local Error \n", "The Error for the 2 step Adams Bashforth is:\n", "\\begin{equation}y_{n+1}=y_n+\\frac{h}{2}[3f(t_{n},w_{n})-f(t_{n-1},w_{n-1})] +\\frac{5h^3}{12}y'''(\\eta),\\end{equation}\n", "where $\\eta \\in [t_{n-1},t_{n+1}]$.\n", "\n", "Rearranging the equations gives \n", "\\begin{equation}\\frac{y_{n+1}-y_n}{h}=\\frac{1}{2}[3f(t_{n},w_{n})-f(t_{n-1},w_{n-1})] +\\frac{5h^2}{12}y'''(\\eta).\\end{equation}\n", "\n", "For our specific initial value problem the error is of the form:\n", "\\begin{equation}\\frac{5h^2}{12}y'''(\\eta)=\\frac{5h^2}{12}2e^{-\\eta} \\leq\\frac{5(0.5)^2}{12} 2\\leq 0.208 \\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "XlS3V2-abVQy", "outputId": "254aa474-540b-4615-b99c-07452ce24e04" }, "source": [ "\n", "d = {'time t_i': t, 'Adams Bashforth, w_i': w,'Exact':y,'Error |w-y|':np.round(np.abs(y-w),5),'LTE':round(2*0.5**2/12*5,5)}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": null, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_iAdams Bashforth, w_iExactError |w-y|LTE
00.01.0000001.0000000.000000.20833
10.50.7130610.7130610.000000.20833
21.00.8032650.7357590.067510.20833
31.51.0040820.9462600.057820.20833
42.01.3268371.2706710.056170.20833
\n", "
" ], "text/plain": [ " time t_i Adams Bashforth, w_i Exact Error |w-y| LTE\n", "0 0.0 1.000000 1.000000 0.00000 0.20833\n", "1 0.5 0.713061 0.713061 0.00000 0.20833\n", "2 1.0 0.803265 0.735759 0.06751 0.20833\n", "3 1.5 1.004082 0.946260 0.05782 0.20833\n", "4 2.0 1.326837 1.270671 0.05617 0.20833" ] }, "metadata": { "tags": [] }, "execution_count": 11 } ] }, { "cell_type": "code", "metadata": { "id": "P_8H3vglbVQ0" }, "source": [ "" ], "execution_count": null, "outputs": [] } ] }